# Import data data<-read.csv("Nestling growth all data FINAL.csv") #str(data) colnames(data)<-c("rowID", "population", "site", "nestID", "age", "date", "time", "nestlingID", "headBill", "culmen", "wingLength", "radius", "tarsus", "tail", "weight", "crop", "flowerScore", "forestCover") data$nestlingID<-as.factor(data$nestlingID) data$site<-as.numeric(data$site) data$nestID<-as.character(data$nestID) data$age<-as.numeric(data$age) # load required library library(nlme) library(lme4) # Modelling approach used below: # For each response variable, calculate a 'global' model. # Then create same for each individual animal # Use the coefs from the global model as starting values # Test whether model parameters vary with site-level canopy cover. # work out which columns to use as response variables responses<-colnames(data)[c(9:15)] n.responses<-length(responses) # remove missing values from these covariates for(i in 1:n.responses) {col<-which(colnames(data)==responses[i]) if(any(is.na(data[, col])==T)==T){data<-data[which(is.na(data[, col])==F),]}} data$population<-as.factor(as.character(data$population)) data$nestlingID<-as.factor(as.character(data$nestlingID)) # remove birds with <=3 obervations birds.all<-xtabs(rep(1, dim(data)[1])~data$nestlingID) birds.excluded<-names(which(birds.all<=3)) for(i in 1:length(birds.excluded)) {data<-data[-which(data$nestlingID==birds.excluded[i]),]} data$nestlingID<-as.factor(as.character(data$nestlingID)) nestlings<-levels(data$nestlingID) n.birds<-length(nestlings) # prepare starting values for global models starting.val.matrix<-matrix(data=NA, nrow=n.responses, ncol=3) colnames(starting.val.matrix)<-c("asym", "xmid", "k") rownames(starting.val.matrix)<-responses starting.val.matrix[, 1]<-c(40, 12, 122.5, 30, 17, 100, 100) starting.val.matrix[, 2]<-c(20, 20, 18.89, 20, 20, 20, 20) starting.val.matrix[, 3]<-starting.val.matrix[, 1]*0.001 starting.val.matrix[3, 3]<-0.136 # import site data #forest<-read.csv("Forest cover per site.csv") # forest$site<-as.character(forest$site) # prepare output objects dataset.list<-list() global.model.list<-list() individual.model.list<-list() predictors<-c("site", "flowerScore", "forestCover", "flowerScore*forestCover") n.predictors<-length(predictors) # prepare dataset to store model reults modelling.output<-as.data.frame(matrix(data=NA, nrow=(6*n.responses* n.predictors), ncol=10)) colnames(modelling.output)<-c( "bird.metric", "response", "predictor", "int.beta", "int.se", "int.P", "predictor.beta", "predictor.se", "predictor.P", "n") for(i in 1:3){modelling.output[, i]<-as.character(modelling.output[, i])} for(i in 4:10){modelling.output[, i]<-as.numeric(modelling.output[, i])} # run a loop that models all response variables for(h in 1: n.responses) { # create global model of wing length starting.vals<-c(asym=starting.val.matrix[h, 1], xmid=starting.val.matrix[h, 2], k=starting.val.matrix[h, 3]) model<-gnls( formula(paste(responses[h], "~ asym/(1 + exp((xmid - age)*k))", sep="")), data=data, start=starting.vals) #summary(model) global.model.list[[h]]<-model # set parameters for individual models starting.vals2<-summary(model)$coefficients #summary(model)$tTable # create output dataframe output<-cbind(data.frame(nestlingID=nestlings), as.data.frame(matrix(data=0, nrow=n.birds, ncol=12))) colnames(output)[2:13]<-c("site", "nest", "flowerScore", "forestCover", "asym", "asym.se", "xmid", "xmid.se", "k", "k.se", "x.min", "x.max") output$nestlingID<-as.character(output$nestlingID) output$site<-as.character(output$site) output$nest<-as.character(output$nest) # run loop to create outputs for individual models for(i in 1:n.birds) { data.thisrun<-data[which(data$nestlingID==nestlings[i]), ] model.thisrun<-try(gnls( formula(paste(responses[h], "~ asym/(1 + exp((xmid - age)*k))", sep="")), data=data.thisrun, start=starting.vals2), silent=T) if(class(model.thisrun)[1]!="try-error"){ output.thisrun<-as.numeric(c( summary(model.thisrun)$tTable[1, 1:2], summary(model.thisrun)$tTable[2, 1:2], summary(model.thisrun)$tTable[3, 1:2], min(data.thisrun$age), max(data.thisrun$age))) output[i, 6:13]<-output.thisrun} # add in details unaltered from 'data' output$flowerScore[i]<-data.thisrun$flowerScore[1] output$forestCover[i]<-data.thisrun$forestCover[1] output$site[i]<-as.character(data.thisrun$population[1]) output$nest[i]<-as.character(data.thisrun$nestID[1]) } # remove errors 1: remove rows where values were not calulated if(any(output$asym==0)==T){output<-output[-which(output$asym==0),]} # remove errors 2: are there any rows where model fit is v. poor? # judge this by SE > (2*beta) remove.rows<-rep(0, dim(output)[1]) test1<-output$asym.se/output$asym if(any(test1>2)==T){remove.rows[which(test1>2)]<-1} test2<-output$xmid.se/output$xmid if(any(test2>2)==T){remove.rows[which(test2>2)]<-1} test3<-output$k.se/output$k if(any(test3>2)==T){remove.rows[which(test3>2)]<-1} output<-output[which(remove.rows==0),] # change relevant vars to factors, save output$site<-as.factor(output$site) output$nest<-as.factor(output$nest) dataset.list[[h]]<-output # create outer loop to include three predictor variables - forest cover, flowering and region for(k in 1: n.predictors) { # create mixed model of all response variables ~ veg cover + (1| nestID/site) for(i in 1:6) { # model for site as predictor (i.e. categorical) if(k==1){model.thisrun<-lme( # if site is a fixed effect, remove as a random effect fixed=formula(paste(colnames(output)[c(6:11)[i]], "~", predictors[k], sep="")), random= ~ 1 | nest, data=output) table.thisrun<-summary(model.thisrun)$tTable row<-which(is.na(modelling.output$bird.metric)==T)[1] modelling.output[row, 4:6]<-as.numeric(table.thisrun[1, c(1,2,5)]) modelling.output[row, 9]<-as.numeric(anova(model.thisrun, test="Chisq")[2, 4]) # replace this model with lme4 to get SEs model.thisrun<-lmer( # if site is a fixed effect, remove as a random effect formula(paste(colnames(output)[c(6:11)[i]], "~", predictors[k], "+ (1|nest)", sep="")), data=output) }else{ # for single, continuous predictors model.thisrun<-lme( fixed=formula(paste(colnames(output)[c(6:11)[i]], "~", predictors[k], sep="")), random= ~ 1 | site/nest, data=output) table.thisrun<-summary(model.thisrun)$tTable row<-which(is.na(modelling.output$bird.metric)==T)[1] if(k==4){ modelling.output[row, 4:9]<-as.numeric( c(table.thisrun[1, c(1,2,5)], table.thisrun[4, c(1,2,5)]))}else{ modelling.output[row, 4:9]<-as.numeric( c(table.thisrun[1, c(1,2,5)], table.thisrun[2, c(1,2,5)]))} # replace this model with lme4 to get SEs model.thisrun<-lmer( # if site is a fixed effect, remove as a random effect formula(paste(colnames(output)[c(6:11)[i]], "~", predictors[k], "+ (1|site/nest)", sep="")), data=output) } # end results for continuous predictors # add in results that don't change with predictor modelling.output$bird.metric[row]<-responses[h] modelling.output$response[row]<-colnames(output)[c(6:11)[i]] modelling.output$predictor[row]<-predictors[k] modelling.output$n[row]<-dim(output)[1] # export each model to a single list individual.model.list[[row]]<-model.thisrun } # end i } # end k print(paste("end h= ", h, sep="")) # check progress } # end h #modelling.output[151:168,] # check outputs signif.rows<-which(modelling.output$predictor.P<=0.05) n.signif<-length(signif.rows) #signif.rows<-signif.rows[c(3,5,7,1,4,6,2)] write.csv(modelling.output[signif.rows, ], # which effects are significant? "curve_ests_significant_ONLY.csv") write.csv(modelling.output, "curve_parameter_significance_FINAL.csv") # export table # plot model predictions from global.model.list quartz(width=10, height=6) par(mar=c(4,4,1,0.5), mfrow=c(2, 4), cex=0.8) # Plot a key x<-seq(-1,1, 0.01) y<-0.9/(1+exp((0-x)*3)) plot(c(0,1)~c(0,1), type="n", ann=F, axes=F, xlim=c(-1, 1), ylim=c(0,1)) box(bty="l") title(xlab="Age", ylab="Response", line=1) lines(x, y) text(-0.45, 0.7, label=expression(paste("y= ", over(a, (1+e^paste("(b-x).k")), sep=""))), cex=1.1) # xmid labels lines(c(0,0), c(0, 0.45), lty=2, col=grey(0.5)) points(0, 0.45, pch=20, bg="black") text(0.02, 0.25, label="b", pos=4) arrows(x0=0.15, x1=0.02, y0=0.2, y1=0.2, code=2, length=0.05, angle=30) arrows(x0=-0.15, x1=-0.02, y0=0.2, y1=0.2, code=2, length=0.05, angle=30) # asymptote labels lines(c(0, 1), rep(0.9, 2), lty=2, col=grey(0.5)) arrows(x0=0.1, x1=0.1, y0=0.98, y1=0.92, code=2, length=0.05, angle=30) arrows(x0=0.1, x1=0.1, y0=0.82, y1=0.88, code=2, length=0.05, angle=30) text(0.1, 0.97, label="a", pos=2) # plot label text(-0.9, 1, label="a)", font=2) # plot model outputs pos1<-rep(2, 7) pos1[6]<-4 x.label<-rep(40, 7) x.label[6]<-0 x<-seq(1, 40, 0.1) for(i in 1:7) { model.thisrun<-global.model.list[[i]] col.thisrun<-which(colnames(data)==responses[i]) coefs.thisrun<-summary(model.thisrun)$coefficients plot(data[, col.thisrun]~data$age, type="n", ann=F, axes=F, ylim=c(0, c(40,12,120,33,16.5,80,100)[i]), xlim=c(0, 40)) axis(1, seq(0, 40, 10)) axis(2, seq(0, c(40,12,120,40,20,80,100)[i], c(10,2,20,10,5,20,20)[i]), las=1) box(bty="l") points(data$age, data[, col.thisrun], pch=19, col=grey(0.5), cex=0.7) y<-coefs.thisrun[1]/ (1+exp((coefs.thisrun[2]-x)* coefs.thisrun[3])) lines(x, y, lwd=2) text(1, c(40,12,120,33,16.5,80,100)[i], labels=paste(letters[(i+1)], ")", sep=""), font=2) title(ylab=c("Head-Bill Length (mm)", "Culmen Length (mm)", "Wing Length (mm)", "Radius Length (mm)", "Tarsus Length (mm)", "Tail Length (mm)", "Weight (grams)")[i], line=2.5) y.max<-c(40,12,120,33,16.5,80,100)[i]*0.25 if(i==6){y.max<-80} y.int<-c(40,12,120,33,16.5,80,100)[i]*0.08 for(j in 1:3) {par<-c("a", "b", "k")[j] val<-round(summary(model.thisrun)$tTable[j, 1], 2) se<-round(summary(model.thisrun)$tTable[j, 2], 2) text(x.label[i], (y.max-(y.int*j)), label=paste(par, " = ", val, " (", se, ")", sep=""), pos=pos1[i])} } par(mfrow=c(1,1)) title(xlab="Age (days)", line=2.5, cex.lab=0.8)