The error message you see says that something is being passed a symbol instead of a number and that it needs a number instead of a symbol. Sometimes it is important to make sure that functions only evaluate when you give them a number and won't evaluate when given them a symbol. You can have functions only evaluate when given a number by using a pattern,?NumericQ. See this piece of documentation here http://support.wolfram.com/kb/12502.

mylik[f0_?NumericQ, p_?NumericQ] := 
  Module[{ans}, 
   ans = -Log[2] - Log[120] - Log[6227020800] - 
     Log[1124000727777607680000] - Log[15511210043330985984000000] + 
     f0 Log[(1 - p)^6] + 25 Log[6 (1 - p)^5 p] + 
     22 Log[15 (1 - p)^4 p^2] + 13 Log[20 (1 - p)^3 p^3] + 
     5 Log[15 (1 - p)^2 p^4] + Log[6 (1 - p) p^5] + 2 Log[p^6] - 
     LogGamma[1 + f0] + LogGamma[69 + f0]];
{maxlog, mles} = NMaximize[{mylik[f0, p], f0 > 0, 0 < p < 1}, {f0, p}];

proflikfun[f0_?NumericQ] := 
  Module[{ans, p}, ans = NMaximize[{mylik[f0, p], 0 < p < 1}, p];
   ans = ans[[1]];
   ans = ans + Quantile[ChiSquareDistribution[1], 0.95]/2 - maxlog;
   ans];

FindRoot uses normally uses things like Newton's method. It differentiates its input which can't be done in this case. It can however use secant like methods which don't use derivatives if you give it two starting points.

FindRoot[proflikfun[f0], {f0, 1, 2}]

NSolve uses a variety of methods but doesn't take a starting point. Instead, you list what variables to solve for. The correct syntax is:

NSolve[proflikfun[f0] == 0, f0]

These functions are numeric solvers. They don't have any way of knowing about the properties of the function proflikfun. So they don't know how many roots it has or even if the function is well behaved. NSolve will try it's best to find all of the roots, but there's guarantee it can. You can always use FindRoot with many different starting values and see if any of them converge. There's no guarantees. Another method might be to rephrase it as a minimization problem and use a global minimizer like NMini

It may be possible to derive a symbolic version of proflikfun. If it's a polynomial or something simple, then there might be a symbolic method to find all of its roots and then Solve might work. Mylik doesn't need to be defined with ?NumericQ or Module. The problem seems to simplify down pretty substantially:

mylik[f0_, p_] := -Log[2] - Log[120] - Log[6227020800] - 
  Log[1124000727777607680000] - Log[15511210043330985984000000] + 
  f0 Log[(1 - p)^6] + 25 Log[6 (1 - p)^5 p] + 
  22 Log[15 (1 - p)^4 p^2] + 13 Log[20 (1 - p)^3 p^3] + 
  5 Log[15 (1 - p)^2 p^4] + Log[6 (1 - p) p^5] + 2 Log[p^6] - 
  LogGamma[1 + f0] + LogGamma[69 + f0]

For a fixed value of f0, when is mylik maximized?

FullSimplify[D[mylik[f0, p], p] == 0]

For a given value of f0, which value of p maximizes mylik?

Solve[FullSimplify[D[mylik[f0, p], p] == 0], p]
{{p -> 145/(6 (68 + f0))}}

So we can get rid of a lot of the numerical algorithms in this code. That might be helpful.