Sum up uncertainty from confidence intervals in R

Question

I'm fitting a linear regression model in R:

head(data)  
   `x1`         `x2`     `y`  
    4.9711     32.70   0.2810632  
    6.756072   31.60   0.3115225  
    5.895213   56.50   0.4853171  
    1.951329   29.80   0.3010127  

EDIT: fit <- glm(y ~ x1 + x2,family = binomial(link = logit), data = data, control= list(epsilon = 0.0001, maxit = 50, trace = F)))

where y is a proportion.

Here's my code:

fit.l.hat=predict.glm(fit,se.fit=TRUE,type="response")    
ci=c(fit.l.hat$fit - 1.96*fit.l.hat$se.fit, fit.l.hat$fit + 1.96*fit.l.hat%se.fit)

For the first fitted value, x1 = 4.97 and x2 = 32.7. So from my above code, I get this output:

Fitted value: 0.35, Lower 95%: 0.01, Upper 95%: 0.68, Standard Error: 0.17

and so on.

So my question is how to actually interpret this for "dummies". I would say "For an input of 4.97 and 32.7, 95 out of 100 samples will conclude a proportion between 0.01 and 0.68." Is that correct?

What other information can I get out of this, or what other "dummy" ways are there to present this information? And more importantly, how do I quantify the error?

Answer 1

There are several issues here.

First, as your response is a proportion, and you are trying to fit a logistic regression (glm with family=Binomial) you need to specify sample size that each of those proportions refer to. See How to do logistic regression in R when outcome is fractional (a ratio of two counts)?. Your procedure only makes sense of underlying those proportions are a number of binary trials, and to know how confident we are in the model we need to know how many trials there were.

If in fact the proportions are not the result of binary trials (eg they are the proportion of some continuous variable, like how full a bottle was) you need a different procedure - the one that springs to mind is a quasi-binomial family for your glm. This will allow an extra parameter for the dispersion of the results.

You have to do one of the above things, or your results are meaningless.

Second, assuming you fix the above issue, in the way you are constructing your confidence interval you are assuming that the predicted values of y are normally distributed on their original scale. If you are prepared to assume this, you might as well have fit an ordinary least squares model (glm with family=Gaussian) and be done with it. The problems with this are several, one of them being that the resulting confidence intervals are quite likely to include values outside of [0, 1]. Addressing this sort of issue is why we would use a different glm in the first place as you have.

A better procedure is to use predict() to create predicted point values and standard errors on the linear predictor scale; convert those into confidence intervals you have ie assuming Normality and hence +/- 1.96 X se; and then transform the upper and lower bounds of those confidence intervals with plogis() into proportions. This will give you non-symmetric confidence intervals (ie the point estimate won't be in the middle of the upper and lower bounds) but this is quite appropriate when the response is a proportion.

Finally, once you have fixed those first to problems, we come to the question you actually asked - how do I interpret this confidence interval? Actually, the answer to that depends somewhat on your answer to the first question. It's possible that what you really want is a prediction interval.

A confidence interval takes into account the uncertainty in the sampling and model fitting process and your interpretation isn't too bad - if you gathered this data and fit the model this way many times, 95% of the time a confidence interval created this way would have the true value of the mean proportion at that point.

A prediction interval goes one step further and also takes into account the uncertainty of taking another sample (of what size?) and what that observed proportion would be. To get that right 95% of the time you need to take into account the randomness of the new process too, not just whether your parameter estimates are correct. As how to do this depends very much on where your proportions come from I can't answer that bit yet.

Answer 2

R fitted a linear model. For the input values x1 and x2, the output values is drawn at random from a distribution with mean 0.35 and SE of 0.17. You are correct in saying there is a 95% chance that the output value for those input values will be between 0.01 and 0.68.