rm(list=ls()) getwd() #setwd("C:/Users/Pedro/Documents/Classes/Methods in Ecology II/Demos/") par(mfrow=c (1,1)) dataforstats <- read.table("dataforstats_final_BB121613.txt", header=T) names (dataforstats) library(glmmTMB) library(bbmle) library(mgcv) library(plot3D) sub1 <-subset(dataforstats,!is.na(upland_elev_m )) summary(sub1) par(mfrow=c (1,2)) hist(sub1$macroct,20,main="Macroinvertebrates",xlab="Number of organisms") hist(sub1$fishct,20,main="Fish",xlab="Number of organisms") y<-as.numeric(sub1$macroct) depth <- sub1$depth depth2 <- depth^2 fish <- log(sub1$fishct+1) fish2 <- fish^2 rancho <- rep(1,length(sub1$ranch)) rancho[sub1$ranch=="bir"] <- 2 rancho[sub1$ranch=="pal"] <- 3 rancho[sub1$ranch=="wil"] <- 4 rancho <- factor(rancho) par(mfrow=c (2,2),mar=c(3,4,1,2)) scatter3D((fish[rancho==1]+1),log(depth[rancho==1]),y[rancho==1],theta = 40, phi = 30, xlab="fish",ylab="depth",zlab="macroinvertebrates",main="ald") scatter3D((fish[rancho==2]+1),log(depth[rancho==2]),y[rancho==2],theta = 40, phi = 30, xlab="fish",ylab="depth",zlab="macroinvertebrates",main="bir") scatter3D((fish[rancho==3]+1),log(depth[rancho==3]),y[rancho==3],theta = 40, phi = 30, xlab="fish",ylab="depth",zlab="macroinvertebrates",main="pal") scatter3D((fish[rancho==4]+1),log(depth[rancho==4]),y[rancho==4],theta = 40, phi = 30, xlab="fish",ylab="depth",zlab="macroinvertebrates",main="wil") M1 <- lm(y ~ depth + depth2 + rancho + fish) summary(M1) M2<-glmmTMB(y~ depth + depth2 + fish + fish2 + rancho, data = sub1, ziformula = ~0,family=nbinom2(link = "log")) summary(M2) #M3 <- gam(y ~ s(depth, fx=FALSE, k=-1, bs="cr")+rancho,family=negbin(c(.01,1))) #summary(M3) #gam.check(M3) M4 <- gam(y ~ s(depth, fx=FALSE, k=-1, bs="cr")+ s(fish,fx=FALSE,k=-1,bs="cr")+rancho,family=nb()) summary(M4) AICctab(M1,M2,M4) ##########################################FISHCT######################################### #macroct title <- levels(sub1$ranch) x<-seq(min(depth),max(depth),1) ################################################# par(mfrow=c (1,4),mar=c(3,4,3,2)) for (k in 1:4) { plot(depth,y+1,type="n",ylim=c(.01,max(sub1$macroct)),xlim=c(-1,max(depth)), main=title[k],xlab= "depth (cm)", ylab=("number of moacroinvertebrates")) ############################################ ord <- order(depth[rancho==k & sub1$fishct ==0]) depth_dat <- depth[rancho==k & sub1$fishct ==0] depth_dat <-depth_dat[ord] macro1 <- y[rancho==k & sub1$fishct ==0] macro1 <- macro1[ord] t <- table(macro1,depth_dat) dep <- unique(depth_dat) mac <- unique(macro1) ord <- order(dep) dep <- dep[ord] orm <- order(mac) mac <- mac[orm] for (i in 1: length(dep)) { for (j in 1: length(mac)) { points(dep[i],(mac[j]), pch=16,cex=log(t[j,i])+1.4,col="blue")}} ############################################ ord <- order(depth[rancho==k & sub1$fishct >0]) depth_dat <- depth[rancho==k & sub1$fishct >0] depth_dat <-depth_dat[ord] mac1 <- y[rancho==k & sub1$fishct >0] mac1 <- mac1[ord] t <- table(mac1,depth_dat) dep <- unique(depth_dat) mac <- unique(mac1) ord <- order(dep) dep <- dep[ord] orm <- order(mac) mac <- mac[orm] for (i in 1: length(dep)) { for (j in 1: length(mac)) { points(dep[i],(mac[j]), pch=16,cex=log(t[j,i])+1.4,col="black")}} ############################################ y2=predict(M2,list(depth=x,depth2=x^2,fish=rep(5,length(x)),fish2=rep(5^2,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response") lines(x,(y2),col="red",lwd=2) y3=predict(M1,list(depth=x,depth2=x^2,fish=rep(5,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response") lines(x,(y3),col="green",lwd=2) #y4=predict(M4,list(depth=x,fish=rep(0,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response") #lines(x,(y4),col="blue",lwd=2) y5=predict(M4,list(depth=x,fish=rep(5,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response") lines(x,(y5),col="black",lwd=2) } ################################################### par(mfrow=c(3,3)) res1 <-residuals(M1) res2 <-residuals(M2) #res3 <-residuals(M3) res4 <-residuals(M4) plot((sub1$macroct)+1,res1, ylim =c(-4,18), log="x",ylab="residuals", xlab="macroinvertebrates") plot((sub1$macroct)+1,res2, ylim =c(-4,18), log="x",ylab="residuals", xlab="macroinvertebrates") #plot((sub1$macroct)+1,res3, ylim =c(-4,18), log="x",ylab="residuals", xlab="macroinvertebrates") plot((sub1$macroct)+1,res4, ylim =c(-4,18), log="x",ylab="residuals", xlab="macroinvertebrates") plot(fish+1,res1, ylim =c(-4,18), log="x",ylab="residuals", xlab="fish") plot(fish+1,res2, ylim =c(-4,18), log="x",ylab="residual", xlab="fish") #plot(fish+1,res3, ylim =c(-4,18), log="x",ylab="residuals", xlab="fish") plot(fish+1,res4, ylim =c(-4,18), log="x",ylab="residuals", xlab="fish") plot(rancho,res1,ylim =c(-3,18),ylab="residuals", xlab="ranches") plot(rancho,res2,ylim =c(-3,18),ylab="residuals", xlab="ranches") #plot(rancho,res3,ylim =c(-3,18),ylab="residuals", xlab="ranches") plot(rancho,res4,ylim =c(-3,18),ylab="residuals", xlab="ranches") par(mfrow=c(1,3)) y_predicted1 <- predict(M1,type="response") y_predicted2 <- predict(M2,type="response") #y_predicted3 <- predict(M3,type="response") y_predicted4 <- predict(M4,type="response") tb1 <-table(sub1$macroct,round(y_predicted1,0)) tb2 <-table(sub1$macroct,round(y_predicted2,0)) #tb3 <-table(sub1$macroct,round(y_predicted3,0)) tb4 <-table(sub1$macroct,round(y_predicted4,0)) fishy <- unique(sub1$macroct) o1 <- order(fishy) fishy <-fishy[o1] pfishy1 <- unique(round(y_predicted1,0)) pfishy2 <- unique(round(y_predicted2,0)) #pfishy3 <- unique(round(y_predicted3,0)) pfishy4 <- unique(round(y_predicted4,0)) o2 <- order(pfishy1) o3 <- order(pfishy2) #o4 <- order(pfishy3) o5 <- order(pfishy4) pfishy1 <- pfishy1[o2] pfishy2 <- pfishy2[o3] #pfishy2 <- pfishy2[o4] pfishy3 <- pfishy2[o5] plot(fishy,fishy,type="n",ylim=c(-1,6),xlab="observed",ylab="predicted", xlim=c(-5,20)) for (i in 1:dim(tb1)[1]) { for (j in 1:dim(tb1)[2]) { if (i==1) v <- tb1[i,j]/30+1 if (i>1) v <-log(tb1[i,j])+0.5 points(fishy[i],pfishy1[j]-1,pch=1,cex=(v),col="black") }} plot(fishy,fishy,type="n",ylim=c(-1,6),xlab="observed",ylab="predicted", xlim=c(-5,20)) for (i in 1:dim(tb2)[1]) { for (j in 1:dim(tb2)[2]) { if (i==1) v <- tb2[i,j]/30+1 if (i>1) v <-log(tb2[i,j])+0.5 points(fishy[i],pfishy2[j]-1,pch=1,cex=(v),col="black") }} #plot(fishy,fishy,type="n",ylim=c(1,6),xlab="observed",ylab="predicted", xlim=c(-5,20)) #for (i in 1:dim(tb3)[1]) { # for (j in 1:dim(tb3)[2]) { # if (i==1) v <- tb3[i,j]/30+1 # if (i>1) v <-log(tb3[i,j])+0.5 # points(fishy[i],pfishy3[j],pch=20,cex=(v),col="black") # }} plot(fishy,fishy,type="n",ylim=c(-1,6),xlab="observed",ylab="predicted", xlim=c(-5,20)) for (i in 1:dim(tb4)[1]) { for (j in 1:dim(tb4)[2]) { if (i==1) v <- tb4[i,j]/30+1 if (i>1) v <-log(tb4[i,j])+0.5 points(fishy[i],pfishy4[j]-1,pch=1,cex=(v),col="black") }} ##########################################FISHCT######################################### #macroct & fish title <- levels(sub1$ranch) x<-seq(min(fish),max(fish),0.1) ################################################# par(mfrow=c (1,4)) for (k in 1:4) { plot(fish,y+1,type="n",ylim=c(-1,max(sub1$macroct)),xlim=c(-1,max(fish)), main=title[k],xlab= "log(fish+1)", ylab=("number of moacroinvertebrates")) ############################################ ord <- order(fish[rancho==k]) fish_dat <- fish[rancho==k] fish_dat <-fish_dat[ord] mac <- y[rancho==k] mac <- mac[ord] t <- table(mac,fish_dat) fisha <- unique(fish_dat) mac <- unique(mac) ord <- order(fisha) fisha <- fisha[ord] orc <- order(mac) mac <- mac[orc] for (i in 1: length(fisha)) { for (j in 1: length(mac)) { points(fisha[i],(mac[j]), pch=16,cex=log(t[j,i])+1.4,col="blue")}} ############################################ y2=predict(M2,list(depth=rep(10,length(x)),depth2=rep(10^2,length(x)),fish=x,fish2 =x^2,rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response") lines(x,(y2),col="red",lwd=2) y3=predict(M1,list(depth=rep(10,length(x)),depth2=rep(10^2,length(x)),fish=x,rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response") lines(x,(y3),col="green",lwd=2) #y4=predict(M4,list(depth=x,fish=rep(0,length(x)),rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response") #lines(x,(y4),col="blue") y5=predict(M4,list(depth=rep(10,length(x)),fish=x,rancho=factor(rep(k,length(x)),levels=levels(rancho))),type="response") lines(x,(y5),col="black",lwd=2) } ################################################### par(mfrow=c(2,2),mar=c(4,4,4,4)) q<-seq(min(depth),41,1) y1<-y2<-y3<-y4 <- array(0,c(40,40)) for (i in 1:length(x)){ for (j in 1:length(q)){ y1[i,j]=predict(M4,list(depth=q[i],fish=x[j],rancho=factor("1",levels=levels(rancho))),type="response") y2[i,j]=predict(M4,list(depth=q[i],fish=x[j],rancho=factor("2",levels=levels(rancho))),type="response") y3[i,j]=predict(M4,list(depth=q[i],fish=x[j],rancho=factor("3",levels=levels(rancho))),type="response") y4[i,j]=predict(M4,list(depth=q[i],fish=x[j],rancho=factor("4",levels=levels(rancho))),type="response") }} image(q,x,y1,xlab="depth",ylab="fish",main="ald",zlim=c(0,4), col=gray((0:10)/10)) image(q,x,y2,xlab="depth",ylab="fish",main="bir",zlim=c(0,4), col=gray((0:10)/10)) image(q,x,y3,xlab="depth",ylab="fish",main="pal",zlim=c(0,4), col=gray((0:10)/10)) image(q,x,y4,xlab="depth",ylab="fish",main="wil",zlim=c(0,4), col=gray((0:10)/10))