###Part I: Preparing the data rm(list=ls()) ##Set your working directory manually or using the following commands## #getwd() #setwd("pathway") #setwd("C:/Users/Pedro/Documents/Classes/Methods in Ecology II/Demos/") ## read the Hypericum data from file orig_data <- read.table("hypericum_data_94_07.txt", header=T) ## change variables to include only reproductive individuals dt <- subset(orig_data, !is.na(ht_init) & rp_init > 0 & year<1997 ) yr <- unique(dt$year) height <- dt$ht_init fruits <- dt$rp_init id <- dt$tag bald <- dt$bald year <- dt$year m1 <- lm(fruits~height) summary(m1) res.M1 <- residuals(m1) par(mfrow=c (1,3)) for (i in 1:3) { plot(height[year==yr[i]],fruits[year==yr[i]], main=yr[i],xlab="height (cm)", ylab="Number of fruits",pch=16,cex=0.55, xlim=c(0,80), ylim= c(0,1000) ) abline(lm(fruits~height), col="red") } par(mfrow=c (1,2)) plot(height,res.M1,pch=16,ylab="residuals", xlab="height", cex=0.5) plot(height[res.M1<500], res.M1[res.M1<500],pch=16, ylab="residuals", xlab="height", cex=0.5, xlim=c(0,80)) m2 <- lm(log(fruits)~log(height)) summary(m2) res.M2 <- residuals(m2) par(mfrow=c (1,3)) for (i in 1:3) { plot(height[year==yr[i]],fruits[year==yr[i]], main=yr[i],xlab="height (cm)", ylab="Number of fruits",pch=16,cex=0.55, xlim=c(0,80), ylim= c(0,1000) ) o <- order(height) pred <- exp(predict(m2,type="response")) lines(height[o],pred[o],pch=16,cex=0.55, col="red") } library(minpack.lm) Power <- nlsLM(fruits~a1*height^a2,start=list(a1=1,a2=1),data=dt, algorithm="LM") summary(Power) res.P <- residuals(Power) AIC(m1,Power) par(mfrow=c (1,3)) for (i in 1:3) { plot(height[year==yr[i]],fruits[year==yr[i]], main=yr[i],xlab="height (cm)", ylab="Number of fruits",pch=16,cex=0.55, xlim=c(0,80), ylim= c(0,1000) ) o <- order(height) pred4 <- predict(Power,type="response") pred <- predict(m2,type="response") lines(height[o],exp(pred[o]),pch=16,cex=0.55, col="red") lines(height[o],pred4[o],pch=16,cex=0.55, col="blue") } par(mfrow=c (1,2)) plot(height,res.M2,pch=16,ylab="residuals", xlab="height", cex=0.5, main="MLS on ln(data)") abline(h=0, col="red") plot(height,res.P,pch=16,ylab="residuals", xlab="height", cex=0.5, main="Likelihhod on data") abline(h=0, col="red") m3 <- lm(log(fruits)~log(height),weights=height) summary(m3) res.M3 <- residuals(m3) par(mfrow=c (1,3)) for (i in 1:3) { plot(height[year==yr[i]],fruits[year==yr[i]], main=yr[i],xlab="height (cm)", ylab="Number of fruits",pch=16,cex=0.55, xlim=c(0,80), ylim= c(0,1000) ) o <- order(height) pred <- exp(predict(m2,type="response")) pred1 <- exp(predict(m3,type="response")) lines(height[o],pred[o],pch=16,cex=0.55, col="red") lines(height[o],pred1[o],pch=16,cex=0.55, col="blue") } AIC(m2,m3) par(mfrow=c (1,3)) hist(res.M1,40) hist(res.M2,40) hist(res.M3,40) r <- table(id,year,bald) s <- rowSums(r[,,13]) s1 <- s[s==3] repeated <- as.numeric(labels(s1)) f<-0 hrep <- rrep <- rep(0,length(repeated)*3) for (j in 1:3){ for (i in 1:length(repeated)){ f <-f+1 dtr <- subset(dt,year==yr[j] & bald ==93) hrep[f] <- dtr$ht_init[dtr$tag==repeated[i]] rrep[f] <- dtr$rp_init[dtr$tag==repeated[i]] } } m4 <- lm(log(rrep)~log(hrep)) summary(m4) res.M4 <- residuals(m4) par(mfrow=c (1,1)) plot(height,fruits, main="Non independent sample",xlab="height (cm)", ylab="Number of fruits",pch=16,cex=0.45, xlim=c(0,80), ylim= c(0,1000),col="gray" ) points(hrep,rrep,col="red",pch=16,cex=0.55) lines(height[o],pred[o],pch=16,cex=0.55, col="black",lwd=2) or <- order(hrep) predr <- exp(predict(m4,type="response")) ################ random iter <- 500 co <-rep(0,iter+1) co[1] <- m4$coeff[2] for(k in 1:iter){ rat <- sample(length(height),length(hrep)) hrep1 <- height[rat] rrep1 <- fruits[rat] m5 <- lm(log(rrep1)~log(hrep1)) summary(m5) co[k+1] <- m5$coeff[2] res.M5 <- residuals(m4) or1 <- order(hrep1) predr1 <- exp(predict(m5,type="response")) lines(hrep1[or1],predr1[or1],pch=16,cex=0.45, col="lightblue",lwd=0.2)} lines(height[o],pred[o],pch=16,cex=2, col="black",lwd=4) lines(hrep[or],predr[or],pch=16,cex=0.35, col="red",lwd=4) hist(co,80,main="",xlab="Coefficient",col="lightblue") abline(v=co[1],col="red") abline(v=m2$coeff[2],col="black",lwd=4)