[R] confidence intervals with glmmPQL

Tue Nov 27 16:23:32 CET 2012

Hi - Thanks for the reply.  I provided answers below.
Sally

Sally_roman <sroman <at> umassd.edu> writes:

>  Hi - I am using R version 2.13.0.  I have run several GLMMs using
> the glmmPQL function to model the proportion of fish caught in one
> net to the total caught in both nets by length.  I started with a
> polynomial regression full model with three length terms: l, l^2,
> and l^3 (l=length).  The length terms and intercept were the fixed
> effects and the random effect was a paired haul (n=18).

> m1<-glmmPQL(fixed=Proportion~1+Length+second+third,random=~1|Pair,
> family=binomial,data=species,verbose=T,niter=2,
> weight=(Experimental+Control))

  Why did you set niter=2?  That seems like a bad idea
(the default is 10; I can imagine increasing it if there
 are warnings that the fit hasn't converged, but I don't
 see why you would decrease it).

I set niter=2 after I ran the models with the default.  All models converged
within 2 iterations.  

> For the majority of the models, I ended up with a constant model with no
> length effect.

This isn't quite clear: did you do some kind of stepwise model
reduction or AIC-based model selection?

I used a step wise backward selection for length terms and terms were
removed if they were not significant.  

> The issue I am having is with the confidence intervals that
> were calculated.  For two models the CIs are not symmetrical around the
> mean
> proportion from the model.  The CIs for the other constant models are
> symmetrical around the mean.  I was wondering if anyone has an idea why
> this
> would be or if anyone has any suggestions.  
> Thanks Sally

  There's not enough detail here to answer.  How did you get your
confidence intervals?  A reproducible example would be helpful.
If m1 is the result of a glmmPQL fit, then confint(m1) gives odd
results (because it calls confint.default, which doesn't really
know what to do with the results of the fit); intervals(m1) might
be more what you're looking for.

#subsample of data
Pair	Species	Length	Experimental	Control	Proportion	second	third
2	Monkfish	23	0	1	0	529	12167
2	Monkfish	27	0	1	0	729	19683
7	Monkfish	27	0	2	0	729	19683
8	Monkfish	27	0	1	0	729	19683
14	Monkfish	28	1	0	1	784	21952
13	Monkfish	29	0	1	0	841	24389
14	Monkfish	29	0	1	0	841	24389
1	Monkfish	30	0	1	0	900	27000
5	Monkfish	30	0	2	0	900	27000
18	Monkfish	30	1	0	1	900	27000
4	Monkfish	31	1	0	1	961	29791
5	Monkfish	31	0	1	0	961	29791
11	Monkfish	31	0	1	0	961	29791
2	Monkfish	32	0	1	0	1024	32768
11	Monkfish	32	0	1	0	1024	32768
12	Monkfish	32	0	1	0	1024	32768
18	Monkfish	32	1	1	0.5	1024	32768
4	Monkfish	34	0	1	0	1156	39304
17	Monkfish	34	0	1	0	1156	39304
18	Monkfish	34	0	1	0	1156	39304
18	Monkfish	35	0	1	0	1225	42875
1	Monkfish	37	0	1	0	1369	50653
2	Monkfish	40	1	0	1	1600	64000
14	Monkfish	40	0	1	0	1600	64000
18	Monkfish	40	1	0	1	1600	64000

This is my model code
#Models
#Select trip
t<-t1
#select species
species<-subset(t,Species=="Monkfish") 
#define polynomials
#2nd order poly
second<-species$Length^2
#3rd order poly
third<-species$Length^3
#combine species,second & third
species<-cbind(species,second,third)

#backward selection
#full model
m1<-glmmPQL(fixed=Proportion~1+Length+second+third,random=~1|Pair,family=binomial,data=species,verbose=T,niter=2,weight=(Experimental+Control))
summary(m1)
#2nd order
m2<-update(m1, . ~. -third)
summary(m2)
#linear
m3<-update(m2, . ~. -second)
summary(m3)
#constant 
m4<-update(m3, . ~. -Length)
summary(m4)

I used some code that was passed to me from a colleague who got it from
someone else.  I tried to go through the code to understand it, but got
stuck on a couple of areas.  I did not hear back from the original person
who wrote the code.   

Here is the code for a constant model - 

#CIs
get.sel.and.conf.band.catch=function(parm,varm,l.min,l.max,n=200){
  lgt<-seq(l.min,l.max,length=n)
  X=sapply(1:length(parm),function(i){lgt^(i-1)})
  reg.line=X%*%parm
  se=apply(X,1,function(x,varm){sqrt(t(x)%*%varm%*%x)},varm)

cbind(lengtht=lgt,min.L50=reg.line-2*se,mean.L50=reg.line,max.L50=reg.line+2*se)
}

#plots model results
plot.fit=function(glmm.fit,limx,Spec.Text){
  varm=glmm.fit$varFix
  coeff=glmm.fit$coeff$fixed
  eta.matrix=get.sel.and.conf.band.catch(coeff,varm,min(limx),max(limx))
  xxx=c(eta.matrix[,1],rev(eta.matrix[,1]))
  yyy=c(ilogit(eta.matrix[,2]),rev(ilogit(eta.matrix[,4])))
  plot(eta.matrix[,1],ilogit(eta.matrix[,3])
	,type="l"
	,xlab="Length (cm)"
	,ylab="Proportion: exp/(exp+control)"
	,ylim=c(0,1), xlim=limx
	,col="black",lwd=2,las=1
	,axes=FALSE
	,family="serif")
  axis(1,labels=T,cex.axis=0.9)
  #axis(3,labels=TRUE,cex.axis=0.9)
  axis(2,labels=TRUE,las=2,cex.axis=0.9)
  box()
  lines(eta.matrix [,1],ilogit(eta.matrix [,2]),col=4)
  lines(eta.matrix[,1],ilogit(eta.matrix[,4]),col=2)
  polygon(xxx,yyy,col="grey",density=NULL,border="black")
  lines(eta.matrix[,1],ilogit(eta.matrix[,3]),col=1,lwd=2)
  abline(h=0.5,lty=2)
  text(sum(limx)/2,0.9,Spec.Text,cex=1)
}

#model plot code
plot.fit(m4,c(min(species$Length),max(species$Length)),paste(unique(species$Species)))

I looked at other code online from Daniel Hocking which is based on your
code and your code from glmm.wikidot.com.  There are differences between
what I received and your code.  

  Since you are fitting a binomial model, you may be looking for
confidence intervals on the linear predictor (logit) scale, in
which case they should always be symmetric around the estimate,
or on the response (probability) scale, in which case they
should always be asymmetric.

  You may want to send questions about mixed models to
the [hidden email] mailing list. 

Thank you - I will send an email to the mixed model mailing list.  

Thanks again

--
View this message in context: http://r.789695.n4.nabble.com/confidence-intervals-with-glmmPQL-tp4649637p4650972.html
Sent from the R help mailing list archive at Nabble.com.