3
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ If your response is a proportion, for logistic regression to make sense you also need to specify what number of cases each line represents. See stats.stackexchange.com/questions/26762/…. Otherwise you should use a quasi-binomial family for your glm. $\endgroup$ – Peter Ellis Jul 20 '13 at 0:51
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you for your descriptive response! Unfortunately, I have tried to tell the people "above" me that I didn't have a whole lot of faith in the fitted model, but I've been told several times, I am here to QUANTIFY the error, NOT reduce it. So, that's unfortunate. HOWEVER, for my own knowledge and learning, I have edited my glm() call a bit by the way of what my colleagues have, but how would I go about specifying the number of trials? In addition, my x1 variable is normally distributed; however, my x2 is a bit skewed. $\endgroup$ – Darla Jul 22 '13 at 12:20
  • $\begingroup$ Also, would this be more appropriate then? glm(y ~ x1 + x2,family = quasi(link = logit), data = data, control= list(epsilon = 0.0001, maxit = 50, trace = F)) $\endgroup$ – Darla Jul 22 '13 at 13:37
  • $\begingroup$ Also, I recieve the same results using predict() as I do using predict.glm() $\endgroup$ – Darla Jul 22 '13 at 14:12
  • $\begingroup$ Your suggested call here is more appropriate than in the original question, if the "proportion" y is the result of something other than yes/no trials at each combination of x1 and x2. If it is is the result of a number of trials at each combination of x1 and x2, someone must have that knowledge, and you make that number a variable in your data data frame and set the weights= argument to point to that variable. Yes, predict() on a glm object just calls predict.glm() - sorry, should have made that clear. $\endgroup$ – Peter Ellis Jul 23 '13 at 6:45
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.