0

Computation of confidence intervals for odds ratios of a glm object is straight forward:

mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
# set up logistic model
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")
# calc odds ratios and confidence intervals
exp(cbind(OR = coef(mylogit), confint(mylogit)))

Dataset and code credit

However, the same approach does not work on glmmPQL objects (with random effects) returned by MASS:glmmPQL().

MWE from glmmPQL() help:

library(nlme) # will be loaded automatically if omitted
mylogit <- summary(glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,
+                 family = binomial, data = bacteria))
exp(cbind(OR = coef(mylogit), confint(mylogit)))

Manually computing the confidence intervals (coefficient + 2* Std. Error) using the coefficients of the summary output also failed for me as I could not find them in the summary object (but they are there when printing the summary to console). Taking them out manually is no longterm approach...

How does this work efficiently for glmmPQL objects?

Community
  • 1
  • 1
pat-s
  • 5,992
  • 1
  • 32
  • 60
  • Can you please include data and/or code that will provide us with a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) ? – Ben Bolker Aug 02 '16 at 15:11
  • A quick look at the `glmmPQL` function tells me that it is estimated using a Penalised Likelihood. And to quote [this vignette](https://cran.r-project.org/web/packages/penalized/vignettes/penalized.pdf) : "standard errors are not very meaningful for strongly biased estimates such as arise from penalized estimation methods". – StatMan Aug 02 '16 at 15:39
  • @Marcel10 Neglecting the usefulness of the standard errors of penalised likelihood estimations, I would like to include them in the final output since I also included the standard errors of the model coefficients (which represent the log odds), i.e. for completeness. – pat-s Aug 02 '16 at 15:56
  • @pat-s I understand, but I am not sure if there is even a consensus in the academic world about how to define standard errors in the case of penalized regressions. A PhD friend of mine told me that this is not the case for LASSO/Ridge regression (I asked him why the `glmnet` package didn't provide standard errors or p-values). I can't say anything for certainty about your method, but I am afraid it is the same in your case (otherwise it would probably be trivial to extract std. errors). So I think (emphasis on think) you can't get them, because they aren't defined (yet..). Anyway, good luck! – StatMan Aug 02 '16 at 19:52
  • @pat-s on the other hand, that quote of my first comment implies that there is definition... Maybe it can serve as a new lead in your search. – StatMan Aug 02 '16 at 19:54

0 Answers0