################################################### ### chunk number 1: utils ################################################### #line 25 "Banta_trondheim.Rnw" ## setwd("~/Documents/Projects/Mixed Models/GLMM") options(width=80, contrasts=c("contr.treatment", "contr.poly"), continue=" ", SweaveHooks=list(fig=function() par(bty="l",las=1))) ### Load the \code{lme4} package (utility functions to unload packages when switching to the other) if (require(cacheSweave)) { setCacheDir("banta_trondheim_cache") } ################################################### ### chunk number 2: libs ################################################### #line 47 "Banta_trondheim.Rnw" library(lme4) library(bbmle) ## for AICtab library(coefplot2) library(reshape) library(gridExtra) library(MCMCglmm) library(ggplot2) source("glmm_funs.R") ## utility functions ################################################### ### chunk number 3: getdata ################################################### #line 80 "Banta_trondheim.Rnw" dat.tf <- read.csv("Banta_TotalFruits.csv") str(dat.tf) ################################################### ### chunk number 4: ################################################### #line 95 "Banta_trondheim.Rnw" dat.tf <- transform(dat.tf, X=factor(X), gen=factor(gen), rack=factor(rack), nutrient=factor(nutrient), amd=factor(amd,levels=c("unclipped","clipped"))) ################################################### ### chunk number 5: oldway eval=FALSE ################################################### ## #line 106 "Banta_trondheim.Rnw" ## dat.tf$X <- factor(dat.tf$X) ## dat.tf$gen <- factor(dat.tf$gen) ## dat.tf$rack <- factor(dat.tf$rack) ## dat.tf$nutrient <- factor(dat.tf$nutrient) ################################################### ### chunk number 6: reps ################################################### #line 133 "Banta_trondheim.Rnw" (reptab <- with(dat.tf, table(popu, gen))) ################################################### ### chunk number 7: latticeplots1 ################################################### #line 167 "Banta_trondheim.Rnw" bwplot(log(1+total.fruits)~interaction(nutrient,amd)|reg,data=dat.tf) dotplot(log(1+total.fruits)~interaction(nutrient,amd)|reg,data=dat.tf) ################################################### ### chunk number 8: ggplots1 ################################################### #line 174 "Banta_trondheim.Rnw" qplot(interaction(nutrient,amd),log(1+total.fruits), data=dat.tf,geom="boxplot")+facet_wrap(~reg,nrow=1) qplot(interaction(nutrient,amd),log(1+total.fruits), data=dat.tf)+facet_wrap(~reg,nrow=1)+stat_sum() ################################################### ### chunk number 9: ################################################### #line 196 "Banta_trondheim.Rnw" ## use within() to make multiple changes within a data frame dat.tf <- within(dat.tf, { gna <- interaction(gen,nutrient,amd) gna <- reorder(gna, total.fruits, mean) }) ################################################### ### chunk number 10: bwplots ################################################### #line 207 "Banta_trondheim.Rnw" print(bwplot(log(total.fruits + 1) ~ gna, data=dat.tf, scales=list(x=list(rot=90)) ## rotate x-axis labels )) ################################################### ### chunk number 11: eval=FALSE ################################################### ## #line 226 "Banta_trondheim.Rnw" ## ggplot(dat.tf,aes(x=gna,y=log(1+total.fruits)))+geom_boxplot()+coord_flip()+ ## theme_bw() ################################################### ### chunk number 12: checklambda ################################################### #line 236 "Banta_trondheim.Rnw" grpVarL <- with(dat.tf, tapply(log(1+total.fruits), list(gna), var) ) summary(grpVarL) ################################################### ### chunk number 13: checklambda ################################################### #line 245 "Banta_trondheim.Rnw" grpMeans <- with(dat.tf, tapply(total.fruits, list(gna), mean)) summary(grpMeans) ################################################### ### chunk number 14: logvar ################################################### #line 273 "Banta_trondheim.Rnw" grpVars <- with(dat.tf, tapply(total.fruits, list(gna), var) ) ## variance of UNTRANSFORMED data ################################################### ### chunk number 15: ################################################### #line 281 "Banta_trondheim.Rnw" ddply(dat.tf,"gna", function(x) with(x,data.frame(mean=mean(total.fruits),var=var(total.fruits)))) recast(subset(dat.tf, select=c(gna,total.fruits)), id.var=1, gna~..., fun.agg=function(x) c(mean=mean(x),var=var(x))) ################################################### ### chunk number 16: ################################################### #line 299 "Banta_trondheim.Rnw" lm1 <- lm(grpVars~grpMeans-1) ## `quasi-Poisson' fit phi.fit <- coef(lm1) lm2 <- lm((grpVars-grpMeans)~I(grpMeans^2)-1) k.fit <- 1/coef(lm2) ## alternatively (more understandable but computationally harder): ## nls1 <- nls(grpVars~grpMeans*(1+grpMeans/k), ## start=list(k=1)) ################################################### ### chunk number 17: logvarplot ################################################### #line 311 "Banta_trondheim.Rnw" plot( grpVars ~ grpMeans, xlab="group means", ylab="group variances" ) abline(c(0,1), lty=2) text(105,500,"Poisson") curve(phi.fit*x, col=2,add=TRUE) ## bquote() is used to substitute numeric values ## in equations with symbols text(110,3900, bquote(paste("QP: ",sigma^2==.(round(phi.fit,1))*mu)), col=2) curve(x*(1+x/k.fit),col=4,add=TRUE) text(104,7200,paste("NB: k=",round(k.fit,1),sep=""),col=4) ## loess fit Lfit <- loess(grpVars~grpMeans) mvec <- 0:120 lines(mvec,predict(Lfit,mvec),col=5) ################################################### ### chunk number 18: ################################################### #line 330 "Banta_trondheim.Rnw" ggplot(data.frame(grpMeans,grpVars), aes(x=grpMeans,y=grpVars))+geom_point()+ geom_smooth(colour="blue",fill="blue")+ geom_smooth(method="lm",formula=y~x-1,colour="red",fill="red")+ stat_function(fun=function(x) x*(1+x/k.fit), colour="purple",fill="purple") ################################################### ### chunk number 19: plot1 ################################################### #line 372 "Banta_trondheim.Rnw" print(stripplot(log(total.fruits+1) ~ amd|popu, dat.tf, groups=nutrient, auto.key=list(lines=TRUE, columns=2), strip=strip.custom(strip.names=c(TRUE,TRUE)), type=c('p','a'), layout=c(5,2,1), scales=list(x=list(rot=90)), main="Panel: Population")) ################################################### ### chunk number 20: plot2 ################################################### #line 394 "Banta_trondheim.Rnw" print(stripplot(log(total.fruits+1) ~ amd|nutrient, dat.tf, groups=gen, strip=strip.custom(strip.names=c(TRUE,TRUE)), type=c('p','a'), ## points and panel-average value -- ## see ?panel.xyplot scales=list(x=list(rot=90)), main="Panel: nutrient, Color: genotype")) ################################################### ### chunk number 21: eval=FALSE ################################################### ## #line 405 "Banta_trondheim.Rnw" ## ggplot(dat.tf,aes(x=amd,y=log(total.fruits+1),colour=nutrient))+ ## geom_point()+ ## ## need to use as.numeric(amd) to get lines ## stat_summary(aes(x=as.numeric(amd)),fun.y=mean,geom="line")+ ## theme_bw()+opts(panel.margin=unit(0,"lines"))+ ## facet_wrap(~popu) ## ggplot(dat.tf,aes(x=amd,y=log(total.fruits+1),colour=gen))+ ## geom_point()+ ## stat_summary(aes(x=as.numeric(amd)),fun.y=mean,geom="line")+ ## theme_bw()+ ## ## label_both adds variable name ('nutrient') to facet labels ## facet_grid(.~nutrient,labeller=label_both) ################################################### ### chunk number 22: GLMfits ################################################### #line 429 "Banta_trondheim.Rnw" glm.lis <- lmList(total.fruits~nutrient*amd|gen,data=dat.tf, family="poisson") ################################################### ### chunk number 23: GLMlis ################################################### #line 450 "Banta_trondheim.Rnw" print(plot(glm.lis,scale=list(x=list(relation="free")))) ################################################### ### chunk number 24: qqmath1 ################################################### #line 466 "Banta_trondheim.Rnw" print(qqmath(glm.lis)) ################################################### ### chunk number 25: ################################################### #line 477 "Banta_trondheim.Rnw" gentab <- with(dat.tf,table(gen)) summary(as.numeric(gentab)) gentab[c("5","6","34")] ################################################### ### chunk number 26: full1 ################################################### #line 486 "Banta_trondheim.Rnw" mp1 <- glmer(total.fruits ~ nutrient*amd + rack + status + (amd*nutrient|popu)+ (amd*nutrient|gen), data=dat.tf, family="poisson") ################################################### ### chunk number 27: overdisp ################################################### #line 504 "Banta_trondheim.Rnw" overdisp_fun(mp1) ################################################### ### chunk number 28: q.overdisp ################################################### #line 511 "Banta_trondheim.Rnw" deviance(mp1) ################################################### ### chunk number 29: full2 ################################################### #line 517 "Banta_trondheim.Rnw" mp2 <- update(mp1,.~.+(1|X)) ################################################### ### chunk number 30: ################################################### #line 528 "Banta_trondheim.Rnw" attr(VarCorr(mp2)$gen,"correlation") ################################################### ### chunk number 31: ################################################### #line 536 "Banta_trondheim.Rnw" printvc(mp2) ################################################### ### chunk number 32: full3 ################################################### #line 545 "Banta_trondheim.Rnw" mp3 <- glmer(total.fruits ~ nutrient*amd + rack + status + (amd+nutrient|popu)+ (amd+nutrient|gen)+(1|X), data=dat.tf, family="poisson") ################################################### ### chunk number 33: ################################################### #line 553 "Banta_trondheim.Rnw" printvc(mp3) ################################################### ### chunk number 34: mpobs ################################################### #line 578 "Banta_trondheim.Rnw" mp_obs <- glmer(total.fruits ~ nutrient*amd + rack + status + (1|X), data=dat.tf, family="poisson") ################################################### ### chunk number 35: ################################################### #line 612 "Banta_trondheim.Rnw" mp4 <- update(mp_obs,.~.+(1|popu)+(1|gen)) ################################################### ### chunk number 36: ################################################### #line 624 "Banta_trondheim.Rnw" overdisp_fun(mp4) ################################################### ### chunk number 37: locscale ################################################### #line 635 "Banta_trondheim.Rnw" locscaleplot(mp4,col=ifelse(dat.tf$total.fruits==0,"blue","black")) ################################################### ### chunk number 38: sdplot ################################################### #line 655 "Banta_trondheim.Rnw" coefplot2(mp4,ptype="vcov",intercept=TRUE,legend=TRUE, legend.x="bottom",xlim=c(0,1.5),cex.pts=1.5) ################################################### ### chunk number 39: reffplot ################################################### #line 693 "Banta_trondheim.Rnw" ## squash margins a bit ... pp <- list(layout.widths=list(left.padding=0, right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) r1 <- ranef(mp4,postVar=TRUE) d1 <- dotplot(r1, par.settings=pp) library(gridExtra) print(grid.arrange(d1$gen,d1$popu,nrow=1)) ################################################### ### chunk number 40: coefplot1 ################################################### #line 715 "Banta_trondheim.Rnw" coefplot2(mp4) ################################################### ### chunk number 41: drop1 ################################################### #line 746 "Banta_trondheim.Rnw" (dd_AIC <- dfun(drop1(mp4))) ################################################### ### chunk number 42: ################################################### #line 749 "Banta_trondheim.Rnw" dd_AIC ################################################### ### chunk number 43: drop1LRT ################################################### #line 752 "Banta_trondheim.Rnw" (dd_LRT <- drop1(mp4,test="Chisq")) ################################################### ### chunk number 44: ################################################### #line 755 "Banta_trondheim.Rnw" dd_LRT ################################################### ### chunk number 45: mp6 ################################################### #line 787 "Banta_trondheim.Rnw" mp5 <- update(mp4, . ~ . - amd:nutrient) anova(mp5,mp4) ################################################### ### chunk number 46: dropamdnut ################################################### #line 793 "Banta_trondheim.Rnw" (dd_AIC2 <- dfun(drop1(mp5))) (dd_LRT2 <- drop1(mp5,test="Chisq")) ################################################### ### chunk number 47: eval=FALSE ################################################### ## #line 813 "Banta_trondheim.Rnw" ## library(doBy) ## pb1 <- PBrefdist(mp4,mp5) ## PBmodcomp(mp4,mp5,ref=pb1) ################################################### ### chunk number 48: ################################################### #line 821 "Banta_trondheim.Rnw" overdisp_fun(mp5) ################################################### ### chunk number 49: reN ################################################### #line 833 "Banta_trondheim.Rnw" reStack <- ldply(ranef(mp5)) print( qqmath( ~`(Intercept)`|.id, data=reStack, scales=list(relation="free"), prepanel = prepanel.qqmathline, panel = function(x, ...) { panel.qqmathline(x, ...) panel.qqmath(x, ...) }, layout=c(3,1))) ################################################### ### chunk number 50: jugglepkgs ################################################### #line 894 "Banta_trondheim.Rnw" detach("package:coefplot2",unload=TRUE) detach("package:lme4",unload=TRUE) ## VarCorr() conflict library(glmmADMB) library(coefplot2) ################################################### ### chunk number 51: eval=FALSE ################################################### ## #line 906 "Banta_trondheim.Rnw" ## gnb2 <- glmmadmb(total.fruits ~ nutrient*amd + ## rack + status, ## random=~(amd*nutrient|gen)+ ## (amd*nutrient|popu), ## data=dat.tf, family="nbinom") ################################################### ### chunk number 52: ################################################### #line 917 "Banta_trondheim.Rnw" load("Banta_glmmADMB_fits.RData") ################################################### ### chunk number 53: ################################################### #line 925 "Banta_trondheim.Rnw" AICtab(fits) ################################################### ### chunk number 54: gnbcoefs ################################################### #line 932 "Banta_trondheim.Rnw" coefplot2(fits[c("gnb1","gnb2")]) ################################################### ### chunk number 55: gnbresid ################################################### #line 942 "Banta_trondheim.Rnw" locscaleplot(fits$gnb1) ################################################### ### chunk number 56: ################################################### #line 951 "Banta_trondheim.Rnw" lapply(VarCorr(fits$gnb1),round,3) ################################################### ### chunk number 57: ################################################### #line 958 "Banta_trondheim.Rnw" lapply(VarCorr(fits$gnb1B),round,3) ################################################### ### chunk number 58: ################################################### #line 972 "Banta_trondheim.Rnw" library(MCMCglmm) ################################################### ### chunk number 59: eval=FALSE ################################################### ## #line 977 "Banta_trondheim.Rnw" ## ## use independent effects (idh); full covariance matrix (us) fails ## m2 <- MCMCglmm(total.fruits ~ nutrient*amd + ## rack + status, ## random=~idh(amd*nutrient):popu+idh(amd*nutrient):gen, ## data=dat.tf, family="poisson",verbose=FALSE) ################################################### ### chunk number 60: ################################################### #line 997 "Banta_trondheim.Rnw" load("Banta_MCMCglmm_fit.RData") summary(m2) ################################################### ### chunk number 61: MCMCglmmtrace ################################################### #line 1010 "Banta_trondheim.Rnw" print(xyplot(as.mcmc(m2$Sol),layout=c(3,3))) ################################################### ### chunk number 62: MCMCglmmvtrace ################################################### #line 1018 "Banta_trondheim.Rnw" print(xyplot(as.mcmc(m2$VCV),layout=c(3,3))) ################################################### ### chunk number 63: MCMCglmmcoefs ################################################### #line 1026 "Banta_trondheim.Rnw" coefplot2(m2) ################################################### ### chunk number 64: MCMCglmmvcoefs ################################################### #line 1030 "Banta_trondheim.Rnw" coefplot2(m2,ptype="vcov",intercept=TRUE) ################################################### ### chunk number 65: ################################################### #line 1036 "Banta_trondheim.Rnw" detach("package:glmmADMB") library(lme4) ################################################### ### chunk number 66: nfit ################################################### #line 1044 "Banta_trondheim.Rnw" lm1 <- lmer(log(1+total.fruits)~nutrient*amd+rack+status+(1|popu) + (1|gen),data=dat.tf) ################################################### ### chunk number 67: nfitdiag ################################################### #line 1056 "Banta_trondheim.Rnw" op <- par(mfrow=c(1,2)) locscaleplot(lm1,col=ifelse(dat.tf$total.fruits==0,"blue", ifelse(dat.tf$total.fruits==1,"purple","black"))) r <- residuals(lm1) qqnorm(r) qqline(r,col=2) par(op) ## restore Original Parameters ################################################### ### chunk number 68: coefcmp ################################################### #line 1067 "Banta_trondheim.Rnw" op <- par(mfrow=c(1,2)) coefplot2(list(lm1,mp4)) coefplot2(list(lm1,mp4),ptype="vcov",intercept=TRUE,xlim=c(0,1.5)) par(op) ################################################### ### chunk number 69: eval=FALSE ################################################### ## #line 1076 "Banta_trondheim.Rnw" ## ## playing around with LMMs some more ## ## ## fit categorical vars as groups rather than 'effects' ## lm3 <- lmer(log(1+total.fruits)~nutrient*amd+rack+status+(1|popu/amd) + (1|gen/amd),data=dat.tf) ## ## var of nested treatment groups == 0 (consistent) ## ## ## try to compute adjusted AIC, compare ## ## sum(log(dT(y)/dy)) ## ## http://www.unc.edu/courses/2006spring/ecol/145/001/docs/lectures/lecture18.htm ## ## T(y)=log(1+y) ## ## T'(y) = 1/(1+y) ## trans_adjust <- -2*sum(log(1/(1+dat.tf$total.fruits))) ## AIC(lm1)+trans_adjust ## AIC(mp4) ## actually a much better fit ... ## AIC(mp4)+2*nrow(dat.tf) ## even if we penalize for observation-specific random effects ... ## ## ## *what* is a worse fit? looking at simulated distributions ## set.seed(101) ## s1 <- trunc(exp(unlist(simulate(lm1)))-1) ## s2 <- unlist(simulate(mp4)) ## mval <- max(c(s1,s2)) ## p1 <- as.numeric(prop.table(table(factor(s1,levels=0:mval)))) ## p1[(max(s1)+1):length(p1)] <- NA ## p2 <- as.numeric(prop.table(table(factor(s2,levels=0:mval)))) ## matplot(0:mval,cbind(p2,p1),type="s",col=1:2,lty=1) ## ## tfun <- function(x) { trunc(exp(x)-1) } ## s1 <- simulate(lm1,1000) ## summary(colMeans(tfun(as.matrix(s1))==0)) ## ## s2 <- simulate(mp4,1000) ## cc2 <- colMeans(as.matrix(s2)==0) ## cc1 <- colMeans(tfun(as.matrix(s1))==0) ## ## mm2 <- apply(s2,2,max) ## mm1 <- apply(tfun(as.matrix(s1)),2,max) ## summary() ## ## ## summary(s1) ## plot(prop.table(unlist(simulate(mp4)))) ## ## about overparameterization ## ## ## trying to fit instead via GLM with log link (get correct AIC?) ## lm2 <- glmer((1+total.fruits)~nutrient*amd+rack+status+(1|popu) + (1|gen), ## family=gaussian(link="log"),data=dat.tf) ## ## ## comparing fitted vs obs for both models ## ff <- rbind(data.frame(m="glmer",f=fitted(mp4),obs=dat.tf$total.fruits), ## data.frame(m="lmer",f=exp(fitted(lm1))-1,obs=dat.tf$total.fruits)) ## ggplot(ff,aes(x=f,y=obs,colour=m))+stat_sum(alpha=0.5)+theme_bw()