Not sure how much this will help but here are a few suggestions.

(1) Put all parameters into the objective function, including the ones being adjusted by the optimization.

(2) Insist that they be explicitly NumberQ to avoid any useless attempt at symbolic preprocessing (by NMaximize or Nintegrate).

(3) Tighten the bounds of integration, say to mu -+5. This seems not to detract from quality of the result and speeds the objective function evaluation.

(4) Maybe roll it down to 400 or so interations in the optimizer. Some experiments indicate that gives a result that does not differ by much from 200 iterations.

I also made the constraints weak since that's going to be done by the code anyway. With these changes, here is what I have.

Clear[mylik]
mylik[data_, kdata_, model_, f0_?NumberQ, mu_?NumberQ, 
   sigma_?NumberQ] := Block[{K0, pj, fj, N0, loglik, above, below},
   K0 = kdata;
   If[model == "LNB", 
    pj = NIntegrate[
         N[PDF[BinomialDistribution[kdata, 1/(1 + Exp[-x])], #]*
           PDF[NormalDistribution[mu, sigma], x]], {x, -Infinity, 
          Infinity}] & /@ Range[0, kdata];];
   fj = Prepend[PadRight[data, kdata], f0];
   N0 = Total[fj];
   above = LogGamma[N0 + 1.];
   below = Total[LogGamma[fj + 1.]];
   loglik = fj.Log[pj];
   loglik = above - below + loglik;
   loglik];

link3 = {679, 531, 379, 272, 198, 143, 99, 67, 46, 32, 22, 14, 9, 5, 
   3, 1};
klink3 = 20;

It takes a little over two minutes. I don't know if that would be viewed as a huge improvement but it seems to be faster than what I started out with.

Timing[
 NMaximize[{mylik[link3, klink3, "LNB", f0, mu, sigma], f0 >= 0, 
   sigma >= 0}, {f0, mu, sigma}, MaxIterations -> 400]]

(* Out[137]= {137.848000, {-47.234921188, {f0 -> 610.928712217, 
   mu -> -2.17379818051, sigma -> 0.987297937367}}} *)

Here is a useful speedup proposed by my colleague Michael Trott. We go back to integrating over the full range for this version.

Clear[mylik2, int]

int[kdata_, Y_, mu_, sigma_, x_] = 
  N[PDF[BinomialDistribution[kdata, 1/(1 + Exp[-x])], Y]*
    PDF[NormalDistribution[mu, sigma], x]];

mylik2[data_, kdata_, model_, f0_?NumberQ, mu_?NumberQ, 
   sigma_?NumberQ] := Block[{K0, pj, fj, N0, loglik, above, below},
   K0 = kdata;
   If[model == "LNB", 
    pj = NIntegrate[
         Evaluate[int[kdata, #, mu, sigma, x]], {x, -Infinity, 
          Infinity}, Method -> "DoubleExponential"] & /@ 
       Range[0, kdata];];
   fj = Prepend[PadRight[data, kdata], f0];
   N0 = Total[fj];
   above = LogGamma[N0 + 1.];
   below = Total[LogGamma[fj + 1.]];
   loglik = fj.Log[pj];
   loglik = above - below + loglik;
   loglik];

link3 = {679, 531, 379, 272, 198, 143, 99, 67, 46, 32, 22, 14, 9, 5, 
   3, 1};
klink3 = 20;

Timing[
 NMaximize[{mylik2[link3, klink3, "LNB", f0, mu, sigma], f0 >= 0, 
   sigma >= 0}, {f0, mu, sigma}, MaxIterations -> 400]]

(* Out[160]= {32.536000, {-47.2342781829, {f0 -> 610.950495085, 
   mu -> -2.17381590688, sigma -> 0.98731594371}}} *)