### non-linear mixed effects models to estimate abs. and ### relative growth rates (AGR and RGR) ### based on code from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2011.00155.x library(nlme) library(ggplot2) library(MASS) library(future.apply) ### ANALYSES ### ## Greenhouse growth ### #setwd("") dat<-read.csv("Common_garden.csv") source("NLME_CI_functions.R") ### RB model ### RB.char.mod.dat<-dat[dat$Species.code=="RB",] ####################### ## LOOP to find best ## ####################### ### ASYMPTOTIC FIT #### ####################### out.df<-data.frame() for(i in levels(RB.char.mod.dat$Soil.type.c)){ temp.dat<-subset(RB.char.mod.dat,Soil.type.c==i) #################################### # Monomolecular fit #################################### tmp.mono <- getInitial(Size ~ SSasymp(Day, Asym, R0, lrc), data = temp.dat) mono.RB.temp<- nlme(Size ~ SSasymp(Day, Asym, R0, lrc), fixed=Asym + R0 + lrc ~ 1 , random=Asym ~ 1|Pot.code.c/Individual, data = temp.dat, start=tmp.mono) AIC.mono<-AIC(mono.RB.temp) fitted <- as.data.frame(mono.RB.temp$fitted)$fixed resid <- resid(mono.RB.temp) mss <- sum((fitted - mean(fitted))^2) rss <- sum(resid^2) R2.mono <- mss/(mss + rss) R2.mono ################################### # Logistic fit ################################### tmp.logis <- getInitial(Size ~ SSlogis(Day, Asym, xmid, scal), data = temp.dat) logis.RB.temp<- nlme(Size ~ SSlogis(Day, Asym, xmid, scal), fixed=Asym + xmid + scal ~ 1 , random=Asym ~ 1|Pot.code.c/Individual, data = temp.dat, start=tmp.logis) AIC.logis<-AIC(logis.RB.temp) fitted <- as.data.frame(logis.RB.temp$fitted)$fixed resid <- resid(logis.RB.temp) mss <- sum((fitted - mean(fitted))^2) rss <- sum(resid^2) R2.logis <- mss/(mss + rss) R2.logis ################################### # Gompertz ################################### tmp.gompertz <- getInitial(Size ~ SSgompertz(Day, Asym, b2, b3), data = temp.dat) gompertz.RB.temp<- nlme(Size ~ SSgompertz(Day, Asym, b2, b3), fixed=Asym + b2 + b3 ~ 1 , random=Asym ~ 1|Pot.code.c/Individual, data = temp.dat, start=tmp.gompertz) gompertz.RB.temp AIC.gompertz<-AIC(gompertz.RB.temp) fitted <- as.data.frame(gompertz.RB.temp$fitted)$fixed resid <- resid(gompertz.RB.temp) mss <- sum((fitted - mean(fitted))^2) rss <- sum(resid^2) R2.gompertz<- mss/(mss + rss) R2.gompertz out<-as.data.frame(cbind(i,R2.mono,R2.logis,R2.gompertz,AIC.mono,AIC.logis,AIC.gompertz)) out.df<-rbind(out.df,out) } out.df ## CONCLUSION: logistic performs best, gompertz just behind ## loop to fit logistic models library(MuMIn) logis.RB=list() AIC.logis=list() R2.logis=list() temp.dat=list() tmp.logis=list() for(i in levels(RB.char.mod.dat$Soil.type.c)){ temp.dat[[i]]<-subset(RB.char.mod.dat,Soil.type.c==i) ################################### # Logistic fit ################################### tmp.logis[[i]] <- getInitial(Size ~ SSlogis(Day, Asym, xmid, scal), data = temp.dat[[i]]) logis.RB[[i]]<- nlme(Size ~ SSlogis(Day, Asym, xmid, scal), fixed=Asym + xmid + scal ~ 1 , random=Asym ~ 1|Pot.code.c/Individual, data = temp.dat[[i]], start=tmp.logis[[i]]) AIC.logis[[i]]<-AIC(logis.RB[[i]]) fitted <- as.data.frame(logis.RB[[i]]$fitted)$fixed resid <- resid(logis.RB[[i]]) mss <- sum((fitted - mean(fitted))^2) rss <- sum(resid^2) R2.logis[[i]] <- mss/(mss + rss) } ######################## ### BOOTSTRAP METHOD ### ######################## sampfun <- function(fitted,data,idvar) { pp <- predict(fitted,levels=1) rr <- residuals(fitted) dd <- data.frame(data,pred=pp,res=rr) ## sample groups with replacement iv <- levels(data[[idvar]]) bsamp1 <- sample(iv,size=length(iv),replace=TRUE) bsamp2 <- lapply(bsamp1, function(x) { ## within groups, sample *residuals* with replacement ddb <- dd[dd[[idvar]]==x,] ## bootstrapped response = pred + bootstrapped residual ddb$Size <- ddb$pred + sample(ddb$res,size=nrow(ddb),replace=TRUE) ## note edited from height to Size return(ddb) }) res <- do.call(rbind,bsamp2) ## collect results return(res) } ### Start loops ### pframe.out=list() output.save<-list() ### LOGISTIC FITS ### ### INCLUDING AGR and RGR estimates set.seed(56565) ## bootstrapping method: https://stats.stackexchange.com/questions/231074/confidence-intervals-on-predictions-for-a-non-linear-mixed-model-nlme for(i in levels(RB.char.mod.dat$Soil.type.c)){ xvals <- with(subset(RB.char.mod.dat,Soil.type.c==i),seq(min(Day),max(Day),length.out=100)) nresamp <- 1100 # not all bootstraps converge, so aim to get > 1000 resamps ## predicted values: useful below # get predicted values from full model pframe <- with(subset(RB.char.mod.dat,Soil.type.c==i),data.frame(Day=xvals)) pframe$Size <- predict(logis.RB[[i]],newdata=pframe,level=0) coef <- fixed.effects(logis.RB[[i]]) params <- transform_param.logis(coef) K = params[1] r = params[2] M0 = params[3] rates = data.frame( times = pframe$Day, M = (M0*K)/(M0+(K-M0)*exp(-r*pframe$Day)), AGR = (r*M0*K*(K-M0)*exp(-r*pframe$Day))/(M0+(K-M0)*exp(-r*pframe$Day))^2) rates$RGRt <- rates$AGR/rates$M rates$RGRm <- r*(1 - rates$M/K) pframe<-cbind(pframe,rates) # bootstrap by sampling and calculate AGR etc for each bootstrapped model fit plan(multiprocess) sampdata<-future_replicate(nresamp,sampfun(logis.RB[[i]],subset(RB.char.mod.dat,Soil.type.c==i),"Individual"),simplify=FALSE) #print(lapply(sampdata,head)) myboot <- function(sampledata) { object <- try(pfun_Size_AGR_RGR(update(logis.RB[[i]],data=sampledata),pframe=pframe),silent=TRUE) if (class(object) != "try-error") return(object) } plan(sequential) output <- lapply(X=sampdata, FUN=myboot) sum(unlist(lapply(output,is.data.frame))) keep<-which(unlist(lapply(output,is.data.frame))) print(i) print("n.success = ") print(sum(unlist(lapply(output,is.data.frame)))) #print(lapply(output,head)) output.dirty<-output output<-output.dirty[keep] output.save[[i]]<-output yvals<-as.data.frame(lapply(output,FUN=function(x){x$Size})) names(yvals)<-c(1:length(output)) CI_boot <- get_CI(yvals,"boot_") AGRvals<-as.data.frame(lapply(output,FUN=function(x){x$AGR})) names(AGRvals)<-c(1:length(output)) CI_boot_AGR <- get_CI(AGRvals,"boot_AGR_") RGRvals<-as.data.frame(lapply(output,FUN=function(x){x$RGRt})) names(RGRvals)<-c(1:length(output)) CI_boot_RGR <- get_CI(RGRvals,"boot_RGR_") pframe.out[[i]] <- data.frame(pframe,CI_boot,CI_boot_AGR,CI_boot_RGR) pframe.out[[i]]$Treatment<-i } pframe.out_df<-do.call("rbind", pframe.out) pframe.out_df$Treatment<-as.factor(pframe.out_df$Treatment)