Can specified nlme model be modified to compare more than 2 populations?

Question

I am attempting to model logistic growth over time for 6 different treatments using nlme. I have been referring to Sofaer et al. 2013 "Advantages of nonlinear mixed models for fitting avian growth curves" for the general approach in R. Key bits of their code are pasted below (see http://www.avianbiology.org/sites/avianbiology.org/files/appendix/jab5719.pdf for original). The key modification I am considering is switching from comparing 2 populations to comparing 6 treatments. I am hoping to check in advance whether the formulation they use to test for differences in fixed population effects in the growth model function [e.g., Asym + Adiff*site] can still work if there were 6 levels of site rather than just 2. [note: The runt effect in their model is not relevant to our needs.]

Comparison of growth trajectories between two populations

Goal is to test for differences in the three logistic growth parameters between two populations. Model included a fixed runt effect on the inflection point that does not differ between populations (called runt)

Create two new dataframes containing data from the Alaskan and Californian population

AKgrowth = subset(NTgrowth, NTgrowth$site == 1) CAgrowth = subset(NTgrowth, NTgrowth$site == 0)

Function including differences between populations and a runt effect on the inflection point:

SiteAKiRUNT = function(Age, site, Younger_1ifknown, Asym, xmid, K, Kdiff, middiff, Adiff, runt){(Asym + Adiff*site)/(1 + exp(((xmid+middiff*site+runt*Younger_1ifknown) - Age)*(K+Kdiff*site)))}

Calculate derivatives:

SiteAKiRUNTDeriv = deriv(body(SiteAKiRUNT)[[2]], namevec = c("Asym", "xmid", "K", "Kdiff", "middiff", "Adiff", "runt"), function.arg= SiteAKiRUNT)

Starting values:

startsiteAKiRUNT = c(Asym = 9, xmid = 3, K = .5, Kdiff=0, middiff=0, Adiff = 0, runt = 0)

Model without random effects:

SiteAKiRUNT_noRE_gnls = gnls(weight_g ~ SiteAKiRUNTDeriv(Age, site, Younger_1ifknown, Asym, xmid, K, Kdiff, middiff, Adiff, runt), data = NTgrowth, start = startsiteAKiRUNT)
summary(SiteAKiRUNT_noRE_gnls)

Syntax for running models mirrors syntax shown above, with updated fixed-effect function. Top-ranked model: Random effects of nest and nestling on the asymptotic mass and the inflection point

SiteAKiRUNT_Ai_NestNestling = nlme(weight_g ~ SiteAKiRUNTDeriv(Age, site, Younger_1ifknown, Asym, xmid, K, Kdiff, middiff, Adiff, runt), fixed = Asym + xmid + K + Kdiff + middiff + Adiff + runt ~ 1, random = Asym + xmid ~ 1 | Nest_ID/Nestling_ID, data = NTgrowth, start = startsiteAKiRUNT)
summary(SiteAKiRUNT_Ai_NestNestling)

I don't believe comparing two treatments at a time is a good solution because:

• calculating differences between all combinations of treatments for each parameter would either take a long time, or it would take a long time to build a function to do so
• comparing a treatment to the next largest and so on is not feasible because the treatments change order depending on parameter, it also limits our ability to compare treatments that are not neighbouring sequentially

Answer 1

Example

library(nlme)

Set seed to reproduce example

set.seed(101)

Simulated data

dat <- expand.grid(time = 1:10, rep = 1:4, trt = letters[1:6])

dat$eu <- with(dat, paste0(trt,"_",rep))

dat$y <- c(replicate(24, SSlogis(1:10, 10, 5, 1.5) + rnorm(10, sd = 0.25)))

Good practice to generate a groupedData for nlme

datG <- groupedData(y ~ time | eu, data = dat)

Visualize

plot(datG)

Fit a logistic to each experimental unit

fitL <- nlsList(SSlogis, datG)

Fit the NLME model

fnm1 <- nlme(fitL, random = pdDiag(Asym + xmid + scal ~ 1))

Extract fixed coefficients

fxf <- fixef(fnm1)

Update the model. fxf[.] is for the intercept for each parameter and the rep(0,5) puts placeholder zeros for the differences

fnm2 <- update(fnm1, fixed = Asym + xmid + scal ~ trt, start = c(fxf[1], rep(0, 5), fxf[2], rep(0, 5), fxf[3], rep(0, 5)))

Test for treatment effects

anova(fnm2)

Residuals look ok

plot(fnm2)