# Non linear regression mixed model

I have a datastet with several fungi strains (20) and I am trying to find out which is the optimal growth temperature for each strain. For that, I put them to growth in plates at different temperatures (22,25,28 and 31ºC).
Below is a reduced dataset with 3 sites:

zz <-(" iso temp diam
Itiquira   22  5.0
Itiquira   22  4.7
Itiquira   22  5.4
Itiquira   25  5.8
Itiquira   25  5.4
Itiquira   25  5.0
Itiquira   28  4.9
Itiquira   28  5.2
Itiquira   28  5.2
Itiquira   31  4.2
Itiquira   31  4.0
Itiquira   31  4.1
Londrina   22  4.5
Londrina   22  5.0
Londrina   22  4.4
Londrina   25  5.0
Londrina   25  5.5
Londrina   25  5.3
Londrina   28  4.6
Londrina   28  4.3
Londrina   28  4.9
Londrina   31  4.4
Londrina   31  4.1
Londrina   31  4.4
Sinop   22  4.5
Sinop   22  5.2
Sinop   22  4.6
Sinop   25  5.7
Sinop   25  5.9
Sinop   25  5.8
Sinop   28  6.0
Sinop   28  5.5
Sinop   28  5.8
Sinop   31  4.5
Sinop   31  4.6
Sinop   31  4.3"
)


I fit a model with a single level factor

diam ~ thy * exp (thq*(temp-thx)² + thc*(temp-thx)³)

# thx: Optimum temperature
# thy: Diameter at optimum
# thq: Curvature
# thc: Skewness

d =  subset(df, iso=="Itiquira")

library(nlme)

de <- groupedData(diam~temp|iso, data=d, order=FALSE)

n0 <- gnls(diam~thy*exp(thq*(temp-thx)^2+thc*(temp-thx)^3),
data=de, start=c(thy=5.5, thq=-0.08, thx=25, thc=-0.01), params=thy+thq+thx+thc~1)
summary(n0)

plot(diam~temp, d)
points(fitted(n1)~I(d$temp+0.25), pch=19) with(as.list(coef(n0)), curve(thy*exp(thq*(x-thx)^2+thc*(x-thx)^3), add=TRUE, col=2))  Is there a way (mixed models?) to fit a single model that accounts for all the fungi strains? As I said, I have particular interest in thx parameter (Optimum temperature) and its confidence interval. Could someone help me with that? Following the reply comment I have the thx parameters, and I am trying to get this output table: # thx lwr upr # Iti 25.13 - - # Lon 24.44 - - # Sinop 26.33 - -  This is a Walmes solution (inbox) that matches with @jaimedash one library(nlme) df <- groupedData(diam ~ temp | iso, data = df, order = FALSE) n0 <- nlsList(diam ~ thy * exp(thq * (temp - thx)^2 + thc * (temp - thx)^3), data = df, start = c(thy = 5.5, thq = -0.01, thx = 25, thc = -0.001)) intervals(n0) #, , thx # lower est. upper # Itiquira 23.43404 25.28318 27.13231 # Londrina 23.81607 24.40439 24.99272 # Sinop 25.28567 26.44975 27.61383  • What about jumping into nlme: n0 <- nlme(diam~thy*exp(thq*(temp-thx)^2+thc*(temp-thx)^3), data=df, fixed=thy+thq+thx+thc~1, random=thy+thq+thx+thc~1|iso, start=c(thy=5.5, thq=-0.08, thx=25, thc=-0.01 Commented Apr 18, 2016 at 17:38 • Great @user20650! Now with: ranef(n0)$thx + fixef(n0)[3] I obtain the estimates of each fungi strain. Do you know how to create the CI for each thx parameter? I edited the question to see the final output table. Commented Apr 18, 2016 at 18:30
• Actually, you can use coef(n0) (and there is also a plot method plot(coef(no)) ). As for the CI's I'm not sure: as (i think) parameters are not necessarily normally distributed - but this really is a question for the stats people. Commented Apr 18, 2016 at 18:48
• Do you want to know the best growth temperature for a 'generic' fungi strain, of which the types listed in iso are just a subpopulation? Or are you trying to find the optimum temperature for these particular strains? (The model output table makes me think the second choice.) Commented Apr 19, 2016 at 0:27
• second choice @jaimedash ... optimum temperature for each particular strains. Edited question. Thanks. Commented Apr 19, 2016 at 12:50

Because you want to predict the optimal temperatures of each strain, treating the strains as fixed makes sense. However, an interaction among the three level factor and the model for the optimum makes for a headache. I don't think there's enough data to fit it.

A random effects model might help. With your full data of 20 groups you may be able to easily fit a large random effects structure as in the other answer (or even larger). But, in that case you don't get much help predicting specific optimum of each strain. (Although you'd estimate the mean parameters of the optimum model more accurately due to partial pooling of data across the three groups.

The simplest approach is to fit a separate model for each strain and predict from that. Here's some code adopted from yours (and also using intervals as in @Walmes answer) that does that and combines the estimates to get the table of estimates and intervals:

dlist =  sapply(levels(df$iso), function(ll) subset(df, iso==ll), simplify = FALSE, USE.NAMES = TRUE) # so we get a named list reslist <- lapply(names(dlist), function(iso_name) { n0 <- gnls(diam ~ thy * exp(thq * (temp - thx ) ^ 2 + thc * (temp - thx) ^ 3), data=dlist[[iso_name]], start=c(thy=5.5, thq=-0.08, thx=25, thc=-0.01)) list(thx = c(iso=iso_name, intervals(n0)$coef["thx", ]), model=n0)
})
data.frame(do.call(rbind, lapply(reslist, function(r) r$thx)))  table:  iso lower est. upper 1 Itiquira 23.076061236415 25.2831326285258 27.4902040206366 2 Londrina 23.7432027069316 24.4043925569263 25.0655824069209 3 Sinop 25.2525659791209 26.4496178512724 27.6466697234239  you can explore the fit of each model using plot(reslist[[1]]$model) or summary(reslist[[1]])

edit oh, as @Walmes pointed out (and the question is currently edited to show) you can do this more cleanly using nlme:nlsList

• Right, but as I said in my question, I have 20 strains and I just gave 3 of them for the reproducible example.. Commented Apr 22, 2016 at 12:37
• I hadn't absorbed that part of the question, clearly Commented Apr 22, 2016 at 19:14

Use intervals() to get confidence intervals of the model parameters. I also think that the mixed model is more appropriate. But a model with random effects in all parameters can be difficult to estimate. With the data you provided, the model is estimated but the correlation between random effects is close to one. I prefer random effects in thy and thx parameters because they are low order parameters and the strains seems to vary most in these parameters.

df <- groupedData(diam ~ temp | iso, data = df, order = FALSE)
str(df)

n0 <- nlme(diam ~ thy * exp(thq * (temp - thx)^2 + thc * (temp - thx)^3),
data = df,
fixed = thy + thq + thx + thc ~ 1,
# random = thy + thq + thx + thc ~ 1 | iso,
random = thy + thx ~ 1 | iso,
start = c(thy = 5.5, thq = -0.08, thx = 25, thc = -0.01))

summary(n0)

plot(augPred(n0, level = 0:1))

intervals(n0, which = "fixed")
ranef(n0)

• Thanks @Walmes.. I appreciate your contribution, specially improving the model structure (I also think there were too much random effects for so many points). plot(augPred(n0, level = 0:1)) worked for you? Not for me... And a final point: do you know how to create the output table at the end of my question? Thanks! Commented Apr 20, 2016 at 10:56
• this is the error when I try to plot the predictive values: Error in sprintf(gettext(fmt, domain = domain), ...) : invalid type of argument[1]: 'symbol' Commented Apr 20, 2016 at 11:02
• nice. augPred is a useful command. @Juanchi Also FWIW I get no errors from augPred after running the commands here Commented Apr 21, 2016 at 0:26
• I agree with your treatment of the ranef, although it's worth noting the between-group correlation in the ranef on thx and thy is still 1. @Juanchi note that the fixef 95% CI for thx is tighter here -- this is what the pooling gets you (possibly at the cost of fitting the data from the particular strain) Commented Apr 21, 2016 at 0:32
• Individual fixed model would be the suitable approach (edited answer detailing it). Thanks @jaimedash and @Walmes! Commented Apr 22, 2016 at 12:34