################################################### ### chunk number 1: ################################################### library(emdbook) library(bbmle) library(MASS) library(nlme) source("chapskel.R") ################################################### ### chunk number 2: ################################################### data(FirDBHFec) firdata = na.omit(FirDBHFec[,c("TOTCONES","DBH","WAVE_NON")]) firdata$TOTCONES = round(firdata$TOTCONES) cones = 1+firdata$TOTCONES lm1 = lm(log(cones)~log(DBH),data=firdata) lm2 = lm(cones~DBH,data=firdata) nbNLL.0 = function(a,b,k) { predcones = a*firdata$DBH^b -sum(dnbinom(firdata$TOTCONES,mu=predcones,size=k,log=TRUE)) } ################################################### ### chunk number 3: ################################################### nbfit.0 = mle2(nbNLL.0,start=list(a=1,b=1,k=1)) cones2 = cones[firdata$DBH>6 & firdata$DBH<8] ndb = fitdistr(cones2,"negative binomial") nln = fitdistr(cones2,"lognormal") ngm = fitdistr(cones2,"gamma") ################################################### ### chunk number 4: ################################################### op=par(mfrow=c(1,2),lwd=2,bty="l",las=1,cex=1.5) plot(firdata$DBH,cones,log="xy", xlab="DBH", ylab="cones+1",col="gray",cex=0.7,axes=FALSE) axis(side=2) axis(side=1,at=seq(5,20,by=5)) box() abline(v=c(6,8),col="darkgray") curve(exp(coef(lm1)[1])*x^coef(lm1)[2],add=TRUE) with(as.list(coef(nbfit.0)),curve(a*(x-1)^b,add=TRUE,lty=2)) curve(coef(lm2)[1]+coef(lm2)[2]*x,add=TRUE,lty=3) par(xpd=NA) legend(10.5,6,lty=1:3, c("power/LN", "power/NB", "lin/normal"), cex=0.9,bty="n") par(xpd=FALSE) plot(density(cones2,adjust=0.8,from=0),ylim=c(0,0.04), main="",xlab="cones+1",lty=2,xlim=c(-20,100), axes=FALSE) axis(side=2) axis(side=1,at=seq(0,100,by=25)) box() text(40,0.04,"6 < DBH < 8") points(cones2,runif(length(cones2),max=0.001),col="gray") with(as.list(ndb$estimate),lines(0:100,dnbinom(0:100,size=size,mu=mu),type="s"), lty=1) with(as.list(nln$estimate),curve(dlnorm(x,meanlog,sdlog),add=TRUE,lty=3)) curve(dnorm(x,mean(cones2),sd(cones2)),add=TRUE,lty=4) abline(v=0,col="darkgray") legend(25,0.038, lty=c(2,1,3,4), c("density","NB","LN","normal"),bty="n") ##with(as.list(ngm$estimate),curve(dgamma(x,shape=shape,rate=rate),add=TRUE,lty=4)) ################################################### ### chunk number 5: eval=FALSE ################################################### ## lm.reg = lm(y~x) ################################################### ### chunk number 6: ################################################### linregfun = function(a,b,sigma) { Y.pred = a+b*x -sum(dnorm(Y,mean=Y.pred,sd=sigma,log=TRUE)) } ################################################### ### chunk number 7: eval=FALSE ################################################### ## mle2(Y~dnorm(mean=a+b*x,sd=sigma),start=...) ################################################### ### chunk number 8: eval=FALSE ################################################### ## lm.poly = lm(y~x+I(x^2)) ################################################### ### chunk number 9: eval=FALSE ################################################### ## lm.mreg = lm(y~x1+x2+x3) ################################################### ### chunk number 10: eval=FALSE ################################################### ## lm.1way = lm(y~f) ################################################### ### chunk number 11: eval=FALSE ################################################### ## aov1fun = function(m1,m2,m3,sigma) { ## Y.pred = c(m1,m2,m3)[f] ## -sum(dnorm(Y,mean=Y.pred,sd=sigma,log=TRUE)) ## } ################################################### ### chunk number 12: eval=FALSE ################################################### ## lm.2way = lm(Y~f1*f2) ################################################### ### chunk number 13: ################################################### aov2fun = function(m11,m12,m21,m22,sigma) { intval = interaction(f1,f2) Y.pred = c(m11,m12,m21,m22)[intval] -sum(dnorm(Y,mean=Y.pred,sd=sigma,log=TRUE)) } ################################################### ### chunk number 14: eval=FALSE ################################################### ## mle2(Y~dnorm(mean=m,sd=sigma),parameters=list(m~f1*f2)) ################################################### ### chunk number 15: eval=FALSE ################################################### ## lm(Y~f*x) ################################################### ### chunk number 16: ################################################### ancovafun = function(i1,i2,slope1,slope2,sigma) { int = c(i1,i2)[f] slope = c(slope1,slope2)[f] Y.pred = int+ slope*x -sum(dnorm(Y,mean=Y.pred,sd=sigma,log=TRUE)) } ################################################### ### chunk number 17: ################################################### data(FirDBHFec) firdata = na.omit(FirDBHFec[,c("TOTCONES","DBH","WAVE_NON")]) firdata$TOTCONES = round(firdata$TOTCONES) cones = 1+firdata$TOTCONES logDBH=log(firdata$DBH) fir.dwi = lm(log(cones)~logDBH*WAVE_NON,data=firdata) fir.dw = lm(log(cones)~logDBH+WAVE_NON,data=firdata) fir.d = lm(log(cones)~logDBH,data=firdata) fir.w = lm(log(cones)~WAVE_NON,data=firdata) fir.0 = lm(log(cones)~1,data=firdata) ################################################### ### chunk number 18: ################################################### op=par(lwd=2,bty="l",las=1,cex=1.5,mar=c(5,4,2,4)+0.1,mgp=c(2.5,1,0)) plot(log(cones)~logDBH,pch=as.numeric(firdata$WAVE_NON), xlab="log(DBH)",ylab="log(cones+1)") abline(a=coef(fir.dw)[1],b=coef(fir.dw)[2]) abline(a=coef(fir.dw)[1]+coef(fir.dw)[3],b=coef(fir.dw)[2],lty=2) par(xpd=NA) legend(2.65,1.5,lty=1:2, pch=1:2,c("nonwave","wave"),bty="n") par(xpd=FALSE) par(op) ################################################### ### chunk number 19: eval=FALSE ################################################### ## n1 = nls(y~a*x^b,start=list(a=1,b=1)) ################################################### ### chunk number 20: eval=FALSE ################################################### ## nlsList(TOTCONES~a*DBH^b|WAVE_NON, ## data=firdata,start=list(a=0.1,b=2.7)) ################################################### ### chunk number 21: ################################################### nlsList(TOTCONES~a*DBH^b|WAVE_NON,data=firdata,start=list(a=0.1,b=2.7)) gnls(TOTCONES~a*DBH^b,data=firdata,start=c(0.1,2.7,2.7), params=list(a~1,b~WAVE_NON)) ################################################### ### chunk number 22: eval=FALSE ################################################### ## gnls(TOTCONES~a*DBH^b,data=firdata,start=c(0.1,2.7,2.7), ## params=list(a~1,b~WAVE_NON)) ################################################### ### chunk number 23: ################################################### data(ReedfrogFuncresp) nls(Killed ~ SSmicmen(Initial, a, b), data = ReedfrogFuncresp) ################################################### ### chunk number 24: ################################################### op=par(lwd=2,bty="l",las=1,cex=1.5,mgp=c(2.5,1,0)) firnls = nls(TOTCONES~a*DBH^b,start=list(a=0.1,b=2.7), data=firdata) plot(cones~DBH,data=firdata,ylab="Cones") with(as.list(coef(firnls)),curve(a*x^b,add=TRUE)) curve(exp(coef(fir.d)[1])*x^(coef(fir.d)[2]),add=TRUE,lty=2) legend("topleft",c("power/normal","power/LN"), lty=1:2,bty="n") par(op) ################################################### ### chunk number 25: ################################################### op=par(lwd=2,bty="l",las=1,cex=1.5,mgp=c(2.5,1,0)) data(ReedfrogFuncresp) glm1 = glm(cbind(Killed,Initial-Killed) ~ Initial, family=binomial,data=ReedfrogFuncresp) glm2 = glm(cbind(Killed,Initial-Killed) ~ Initial, family=binomial(link="log"),data=ReedfrogFuncresp) with(ReedfrogFuncresp,plot(Killed/Initial~Initial,xlab="Initial density", ylab="Fraction killed")) ivec = 1:100 lines(ivec,predict(glm1,newdata=list(Initial=ivec),type="response")) lines(ivec,predict(glm2,newdata=list(Initial=ivec),type="response"),lty=2) legend("topright",c("logistic regression","log-binomial model"), lty=1:2,bty="n") par(op) ################################################### ### chunk number 26: ################################################### glm(cbind(Killed,Initial-Killed)~Initial, data=ReedfrogFuncresp,family="binomial") glm(cbind(Killed,Initial-Killed)~Initial, data=ReedfrogFuncresp,family="binomial") ################################################### ### chunk number 27: eval=FALSE ################################################### ## glm1 = glm(y~x,family="poisson") ################################################### ### chunk number 28: ################################################### poisregfun = function(a,b) { Y.pred = exp(a+b*x) -sum(dpois(y,lambda=Y.pred,log=TRUE)) } ################################################### ### chunk number 29: eval=FALSE ################################################### ## glm2 = glm(cbind(y,N-y)~x,family="binomial") ################################################### ### chunk number 30: ################################################### logistregfun = function(a,b) { p.pred = exp(a+b*x)/(1+exp(a+b*x)) -sum(dbinom(y,size=N,prob=p.pred,log=TRUE)) } ################################################### ### chunk number 31: eval=FALSE ################################################### ## glm3 = glm(cbind(y,N-y)~x,family=binomial(link="log")) ################################################### ### chunk number 32: ################################################### logregfun = function(a,b) { p.pred = exp(a+b*x) -sum(dbinom(y,size=N,prob=p.pred,log=TRUE)) } ################################################### ### chunk number 33: eval=FALSE ################################################### ## glm4 = glm(cbind(Killed,Initial-Killed) ~ Initial, ## family=quasi(link="inverse",variance="mu(1-mu)"), ## data=ReedfrogFuncresp) ## x = 10:20 ## p = 1/(1.2+5.5*x) ## y = rbinom(length(x),size=x,prob=p) ## testglm = glm(cbind(y,x-y)~x, ## family=quasi(link="inverse",variance="mu(1-mu)"), ## start=c(1.9,-0.4)) ## testglm = glm(cbind(y,x-y)~x, ## family=quasi(link="inverse",variance="constant"), ## start=c(1.2,5.5)) ################################################### ### chunk number 34: ################################################### a = round(coef(glm1)[1],3) b = round(coef(glm1)[2],4) ################################################### ### chunk number 35: ################################################### glm(TOTCONES~log(DBH),data=firdata,family="quasipoisson") ################################################### ### chunk number 36: ################################################### glm.nb(TOTCONES~log(DBH),data=firdata) ################################################### ### chunk number 37: ################################################### logcones = log(firdata$TOTCONES+1) lm.0 = lm(logcones~1,data=firdata) lm.d = lm(logcones~log(DBH),data=firdata) lm.w = lm(logcones~WAVE_NON,data=firdata) lm.dw = lm(logcones~log(DBH)+WAVE_NON,data=firdata) lm.dwi = lm(logcones~log(DBH)*WAVE_NON,data=firdata) ################################################### ### chunk number 38: ################################################### anova(lm.0,lm.d,lm.dw,lm.dwi) AIC(lm.0,lm.d,lm.w,lm.dw,lm.dwi) ################################################### ### chunk number 39: ################################################### ancovafun = function(i1,i2,slope1,slope2,sigma) { int = c(i1,i2)[WAVE_NON] slope = c(slope1,slope2)[WAVE_NON] Y.pred = int+ slope*log(DBH) -sum(dnorm(logcones,mean=Y.pred,sd=sigma,log=TRUE)) } m1 = mle2(ancovafun,start=list(i1=-2,i2=-2,slope1=2.5,slope2=2.5,sigma=1), data=firdata) AIC(m1) ################################################### ### chunk number 40: ################################################### coef(lm.dwi) coef(m1) ################################################### ### chunk number 41: eval=FALSE ################################################### ## mle2(log(TOTCONES+1)~dnorm(logDBH*WAVE_NON),data=firdata,start=...)