Hi - I am using R version 2.13.0. I have run several GLMMs using the glmmPQL function in the MASS package to model the proportion of fish caught in one net to the total caught in both nets by length for 5 different species. 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). 

For the majority of the models, I ended up with a constant model with no length effect. 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. 

I provided a small sample of data for one species along with the CI, model, and plotting codes. I used some code that was passed to me from a colleague who got it from someone else. I have gone 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. 

#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) 


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 and Ben Bolker at glmm.wikidot.com. There are differences between what I received and the examples from Daniel and Ben. 

Any help would be appreciated. 
Thanks Sally 

-- 

Sally Roman 
Fisheries Biologist 
University of Massachusetts Dartmouth 
School for Marine Science and Technology 
200 Mill Road Suite 325 
Fairhaven, MA 02917 
phone 508-910-6379 

	[[alternative HTML version deleted]]

