I'm trying to fit a generalized mixed effects model to a data set where each subject has paired categorical responses y (so I'm trying to use a binomial logit link). There are about 183 observations and one explanatory factor x. I'm trying to fit something like: (lmer(y~x+(1|subject))) I also tried fitting the same type of model using glmmPQL from MASS. In both cases, I get a t-statistic that is huge (in the thousands) and a tiny p-value. (Just for comparison, if I use lrm ignoring the clustering I get a t-statistic around 3 or so and what appears to be a reasonable estimated coefficient which is very close to the estimated coefficient using just one observation from each subject. Most of the subjects have two responses and in almost all cases the responses are identical although the explantory factor values are not always identical for each subject. If I use geeglm from geepack, I get reasonable estimates close to the naive model results. I also tried using the SAS glimmix macro to fit a generalized mixed model and the routine does not converge. Why does geeglm appear to work but not lmer and glmmPQL? Is this likely to be due to my particular data set? Rick B.
GEE models are very different from the subject-specific models fitted by glmmPQL: see the comparison in MASS. You are testing quite different hypotheses. You appear to be assuming that some things `appear not to work' because they do not give the same results as a different test. Suppose x were logical, and that almost all subjects did better with x=TRUE. Then what you are saying is that most subjects have response 0 or 1 irrespective of x, but suppose that when they differ, it was (almost) always x=TRUE that gave 1. Then the subject-specific model will have a large positive coefficient for x. (It is possible that if the pattern is the same for all individuals the MLE is infinite, called complete separation.) OTOH, the GEE model applies to the population, and in the population x=TRUE makes rather little difference to the mean response. GEE models attenuate subject-specific effects, and can do so dramatically. On Wed, 30 Nov 2005, Rick Bilonick wrote:> I'm trying to fit a generalized mixed effects model to a data set where > each subject has paired categorical responses y (so I'm trying to use a > binomial logit link). There are about 183 observations and one > explanatory factor x. I'm trying to fit something like: > > (lmer(y~x+(1|subject))) > > I also tried fitting the same type of model using glmmPQL from MASS. In > both cases, I get a t-statistic that is huge (in the thousands) and a > tiny p-value. (Just for comparison, if I use lrm ignoring the clustering > I get a t-statistic around 3 or so and what appears to be a reasonable > estimated coefficient which is very close to the estimated coefficient > using just one observation from each subject. > > Most of the subjects have two responses and in almost all cases the > responses are identical although the explantory factor values are not > always identical for each subject. > > If I use geeglm from geepack, I get reasonable estimates close to the > naive model results. > > I also tried using the SAS glimmix macro to fit a generalized mixed > model and the routine does not converge. > > Why does geeglm appear to work but not lmer and glmmPQL? Is this likely > to be due to my particular data set? > > Rick B. > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
>>>>> "Rick" == Rick Bilonick <rab45 at pitt.edu> >>>>> on Wed, 30 Nov 2005 23:11:07 -0500 writes:Rick> I'm trying to fit a generalized mixed effects model to a data set where Rick> each subject has paired categorical responses y (so I'm trying to use a Rick> binomial logit link). There are about 183 observations and one Rick> explanatory factor x. I'm trying to fit something like: Rick> (lmer(y~x+(1|subject))) If you want binomial you have to give ' family = binomial ' to lmer ! Further, always using 'data = ..' is generally recommended practice, and I'd rather assign the result of lmer(.) than just print it, i.e. (if you want to print too): (model1 <- lmer(y ~ x + (1|subject), data = <your DFrame>, family = binomial)) and in your case y should be the 2-column matrix cbind(successes, failures) Martin Maechler, ETH Zurich Rick> I also tried fitting the same type of model using glmmPQL from MASS. In Rick> both cases, I get a t-statistic that is huge (in the thousands) and a Rick> tiny p-value. (Just for comparison, if I use lrm ignoring the clustering Rick> I get a t-statistic around 3 or so and what appears to be a reasonable Rick> estimated coefficient which is very close to the estimated coefficient Rick> using just one observation from each subject. Rick> Most of the subjects have two responses and in almost all cases the Rick> responses are identical although the explantory factor values are not Rick> always identical for each subject. Rick> If I use geeglm from geepack, I get reasonable estimates close to the Rick> naive model results. Rick> I also tried using the SAS glimmix macro to fit a generalized mixed Rick> model and the routine does not converge. Rick> Why does geeglm appear to work but not lmer and glmmPQL? Is this likely Rick> to be due to my particular data set? Rick> Rick B. Rick> ______________________________________________ Rick> R-help at stat.math.ethz.ch mailing list Rick> https://stat.ethz.ch/mailman/listinfo/r-help Rick> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
On Thu, 2005-12-01 at 10:13 +0000, Prof Brian Ripley wrote:> On Thu, 1 Dec 2005, Martin Maechler wrote: > > >>>>>> "Rick" == Rick Bilonick <rab45 at pitt.edu> > >>>>>> on Wed, 30 Nov 2005 23:11:07 -0500 writes: > > > > Rick> I'm trying to fit a generalized mixed effects model to a data set where > > Rick> each subject has paired categorical responses y (so I'm trying to use a > > Rick> binomial logit link). There are about 183 observations and one > > Rick> explanatory factor x. I'm trying to fit something like: > > > > Rick> (lmer(y~x+(1|subject))) > > > > If you want binomial you have to give ' family = binomial ' to lmer ! > > I noticed that, but he did say `something like'. You need to specify the > family for gee and glmmPQL too. > > I think the moral is to give the exact code you use. > > > Further, always using 'data = ..' is generally recommended practice, > > and I'd rather assign the result of lmer(.) than just print it, > > i.e. (if you want to print too): > > > > (model1 <- lmer(y ~ x + (1|subject), data = <your DFrame>, family = binomial)) > > > > and in your case y should be the 2-column matrix > > cbind(successes, failures) > > Not necessarily. If these are binary responses, you can just give y as > shown.Sorry, here is the more complete information: (lmer(y~x+(1|subject),data=mydata,family=binomial)) y consists of zeroes and ones. At the time I was able to post I was working from memory unfortunately. I did use "family=binomial" for all the models. I get the same results whether I assign the results or not. I was just trying to give the basic syntax for the model. Sorry for any confusion. As I said, I think it has to do with the fact that the responses are so highly correlated. The same data fails to converge when using SAS and the glimmix macro (I don't yet have accesss to the new "proc glimmix".) I also made up some artificial data sets and whenever the paired responses were identical the same problem appeared. Unfortunately I can't share the data sets. Do I need to specify the correlation structure explicitly? I thought my data set was similar to others that used the same type of model and functions successfully. Rick B.