Hi Jochen,

When applying confint to a model obtained from glm.nb from the MASS package, it prints "Waiting for profiling to be done..." before any result, which is consistent with base R confint methods for non-linear fits (including fits by glm). Are you certain that the results are not profiled?

The code source of confint.glm indeed calls profile.glm, whose documentation page says it profiles the likelihood.

If you have numerical difference, however, it may be for another reason (according to the doc): the chi-square is asymptotic; the documentation suggests a finite-size correction may be used, since it speaks of a T statistic (but that's a personnal guess only).

Dear Emmanuel,

yes, I am sure that the R function calculates the limits based on the likelihood (using the MLE of the nuisance parameter), not on the profile-likelihood. I can send you a script that demonstrates that.

I think that "profiling" here means to actually find the roots of the logLik - qchisq. But Iam not sure. Maybe to stress the point that the limits are really taken from the logL and not from some quadratic approximation (like in Wald tests, for instance). At least this would fit to the output.

Yet, odds are high that the misconception/error is on my side...

I'm interested by the script... Would help to understand what's going on!

I will prepare it tomorrow (cross check, comment etc). CU

A bit of reflexion, by analogy with (easier) linear model with two parameters - the fact that one is "nuisance" is just a view of mind, since anyway one is thinking around the total likelyhood maximum (profiling being just a trick to find more easily this maximum).

The surface is then a paraboloid, and slices parallel to the parameters plane are ellipses. The ellipses main axis are not colinear with the axes (each axis being for one of the two parameters), because of the correlation between the two parameter estimators.

However, the univariate confidence interval for a parameter is parallel to the axis, and goes through the minimum (it is given by the diagonal term of the matrix defining the ellipsis, in the original basis, not in the rotated basis where the matrix would be diagonal). It is indeed narrower than the one that would be obtained by projecting the ellipsis ends on the interval, because it does not take into account the fact that some parameter couples are not very probable, but that's ok - projecting the end would give a confidence interval with much more confidence than the nominal level. Another way to think this is to say that confidence intervals are "marginal".

Now, if you think to your question in a more general case: if you keep the nuisance parameter constant, you do the analog of using a confidence interval parallel to the parameter axis, that is a marginal one. If you change the nuisance parameter value so has to stay on local minima, then you go along the ellipsis main axis (with the smallest angle with the parameter axis), hence you will reach the end of the ellipsis and project this on the parameter exis, so this should give a wider confidence interval, but with more coverage than required. And not a marginal one, but a somehow "conditional on the fact that the other parameters remains minimum" one.

Hope this is clear, correct (that's just my intuition) and answer your question...

I just prepared the R script with an example for the negative binomial.

Regarding your last answer: I understand this. What puzzles me:

The glms model a mean function (the linear predictor is linked to the expected value), and inferences are about the means. The distributions are thus parametrized to have a "mean" parameter and some of them have a further parameter, that is - to me - in fact a a nuisance parameter (relevant is only the parameter linked to the linear predictor). I am not sure if I am correct here, but this is the way I understand it.

In the Gaussian case, this is done expicitly (via the quantile from the F-distribution; the limits refer to the log-likelihood (not the profile log-Likelihood) with sigma² fixed at the sample variance (not the MLE)). After all, sigma² is treated here as a nuisance parameter. I wonder why this should be different for other glms.

I checked the coverage probability of confint() for the negative binomial model. Depending on the "true" parameter values and the sample size, this is considerably lower than 95% (can be ~80%). The coverage probability based on the profile likelihood is also lower than 95%, but it is not as bad as that of confint().

Thanks for the script, I will check this.

I think the term "nuisance" is just a matter of interpretation, or eventually as a trick to help finding the minimum, but as little to do with the way to compute confidence intervals. After all, in the Gaussian case, the model is exactly the same if you are indeed interested in the variance! And when you have more than one parameter of interest, you still apply the same method to obtain confidence intervals.

Beside, in the Gaussian context, the fact that the "sample" variance (more exactly, the unbiased estimator of the variance) is used instead of the MLE is basically to have more simple formula (avoid the n/(n-k) correction term in all formula), but it is because obtained formula based on the Student's law or the Fisher-Snedecor law are *exact* (which cannot be done for other models) and not asymptotic. Note that there are some book that present the formula using the MLE estimator of the sigma. More difficult to read, in my opinion, but that's another problem.

I think you low coverage by simulation could be precisely because of the finite sample size, whereas chi-square based formula assume an infinite sample size, and the effect of this may change on the convergence speed which depends on the true parameter values. If so, for a given parameter value and assuming you're doing enough simulations to have precise coverage estimates, you should notice that the coverage increases with sample size.

For instance, if you're using a Gaussian quantile instead of a Student one to build a confidence interval, the coverage will be 92 % for 10 degrees of freedom, and 81 % for 2 degrees of freedom.

However, I don't know if using "your" profiled likelihood is a way to correct for finite sample size or may be misleading...

Eventually, would you have the R script as R or Rmd file, and not pdf, to help making tests? Seeing the likelihood profile, with n = 6, I'm pretty sure the problems mainly come from the asymptotic assumption which is not fullfilled (contours are not at all ellipses), hence real coverage could be anything.

Seeing the contour, I would have expected the "fixed size" to be wider, since the contour are very extended along the X-axis, whereas going "in diagonal" (changing) size seems to intersect the contours at x values closest from the maximum. I guess this is because increasing mu is associated with a slight decrease of size, could you make it visible by 1) plotting in red the contour to be intersected and 2) plot the profiled likelihood "minimum mu for a given size" curve on this contour plot? (should be easier for you than from the PDF...).

Thanks for your input. Very much apprechiated!

The low coverage is certainly due to the limited sample size and the fact that Wilks theorem is only asymptotically true. In small sample sizes, the estimate of the nuisance parameter is very uncertain, and this additional uncertainty should be considered in tests / confidence intervals. And I thought using the profile logLik (rather than the logLik evaluated at the MLE of the nuisance parameter) is a way to acknowledge this extra-uncertainty (at least to some extend). Is this wrong? If it is wrong, then why is the profile-Likelihood used at all? And if it is right, I still don't understand why confint() is not using it.

I now attached the R script file (simple script, no markdown).

Hi Jochen, Thanks for the script. I've explored a little bit, going from the glm.nb to an usual lm because the glm.nb has very strong convergence issues for small sample sizes, so that makes simulations very difficult to interpret. I think that there's two things going on, and generalization may be difficult. First, the finite sample sizes makes the classic "horizontal" intervals much too narrow, and the method of the profiled likelihood increases the width of these intervals, hence leading to better coverage. But the correction is far from enough, if you compare them to the exact confidence intervals obtained from the Student law, so that basically does not correct the problem, it seems to be just a side effect (?). Second, in both the N(µ, sigma) and the NB(µ, size) cases, it appears that the two parameters are somehow uncorrelated (the ellipsis are parallel to the axes), so that the correlations effects discussed previously, that were the reason to prefer the "constant other parameters" over the "profiled other parameters" (they being either of interest or nuisance), in practice do not apply in these cases. The question is, is it always true that nuisance parameters will be « orthogonal » to parameters of interest (here, only a case with intercept is explored) or is it a special case here? And when (if possible) it's not true, what has the strongest effect: the overcoverage due to correlations between parameters or the undercoverage due to finite sample size? Of course, when sample size increases, the two intervals are closest and closest, because the shape of the LL surface comes closer to an ellipsis (which would suggest that in such cases, if parameters are correlated, the overcoverage effect would win for intermediate sample sizes). All of these clearly needs either a mathematical proof or further exploration... I attach the scripts in case you want to play... I slightly enriched the displays to have a better understanding of what's going on.

Wow, thanks for that work. It will be quite some work to play with that. It's an interesting question if the two parametes (µ and a scale parameter) are neccesarily uncorrelated. I don't know. The open questions are a bit over my abilities...

However, do you have any knowledge about how confidence intervals based on the likelihood should be calculated? I still think that Pawitan uses a profile-aproach and R uses the likelihood with the "nuisance" parameter set to its MLE. Am I correct here, and if so: what are the reasons for that discrepancy (in this case the coverage of R's CIs is worse - but why should R calculate CIs with clearly suboptimal coverages, at least for small samples)?

I am sure that R is doing something very clever and I don't get it...

Hi Jochen,

Not sure for Pawitan's computation: the text is not quite clear, because it describes the process to *find the minimum* [which makes full sense], but it does not says after that what to do exactly to obtain confidence intervals and so on once the minimum was found. I'll nevertheless call the method you suggest as « Pawitan-CI » for simplicity below.

After a night for clarifying my ideas, my position is now as follows:

- confidence intervals for a given parameter should be computed at all other parameters fixed, because they are « marginals » hence parallel to the axes;

- the real problem is to have the correct contour that must be crossed to have the expected coverage;

- intervals obtained by « Pawitan-profile » are projection of the « ellipsis » main axis on the parameter-of-interest axis, hence are not confidence intervals, or confidence intervals of unknown coverage because it depends on the amount of correlation between the parameters (the ellipsis axis angle with the parameter axis). In fact, they may be could be seen as « main axis » of the *joint confidence region* of the parameters. However, what is certain is that they will be wider than the marginal confidence intervals obtained from the same contour, so they will ve over-confident;

- in the Gaussian case (and, apparently, the binomial negative law) the main axis of the « ellipsis » are parallel to the axis [considering µ and sigma], hence confidence interval are the main axis of the ellipsis that corresponds to the correct quantile. Consequently, « Pawitan-profile » CI and « R-profile CI » should be the same.

For the Gaussian model, it is a proved fact that arithmetic mean and estimator of the variance are independent, hence orthogonals, so this observation was in fact expected. As far as I know, the Gaussian distribution is the only continuous one for which this is true, so this may not be a general result (but the negative binomial is not continuous, so I don't know if the results hold or not). Exploring the estimation of parameters of another law could be quite interesting, by the way, to check this.

Now, in finite samples, two things occur:

1) the chisq contour is too small because it neglects the "nuisance parameters" uncertainty — so the R-profile CI will be undercovered.

2) the shape of the likelihood is « distorted » from an ellipsis, which means that "minimal lines" are not any more straight lines (the ellipsis axis) but curved. Hence, R-profile CI and Pawitan-profile CI are not identical, even if they "should" be.

3) It's not clear why this distortion should always translate into a wider interval, but it happens to be the case in both the lm and the glm.nb cases.

4) It's not clear why this distortion should always improve the coverage of the CI, and why it could not in fact over-correct. For the lm case, it seems to under-correct; for the glm.nb case, the difference between the two is quickly practically irrelevant...

In fact, with increasing sample size, the distortion of the shape may reduce quicker than the convergence rate of the quantile toward the chi-square quantile (that seems to be the case for the Gaussian case), hence the difference between the two intervals becomes irrelevant but both still give under-covering CI. Or, conversely, if the contour shape converges to an ellipsis more slowly than the quantile (which seems less probable, but who knows...), the Pawitan-profile CI will overcorrect for "moderate" sample sizes.

Adding that if nuisance parameters is not orthogonal to the other parameters, it's not clear at all in which direction it will go, I guess the reason to stick to the R-profile CI is that it is logically sounded, despite we know it is undercovered, whereas one does not really know what we are doing when trying to find the Pawitan-profile CI.

In fact, focusing on using the good contour should be a better option than trying to distort the « CI » to make it longer. That is what does the bootstrap, if I am correct.

Of course, in cases the "distortion" approach works, it's a good thing, but I guess unless simulations are done, it's difficult to be certain for all possible models that it is an improvement, which is problematic for a generic function.

But that's just an opinion...