5
$\begingroup$

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"
)
df <- read.table(text=zz, header = TRUE)

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))

enter image description here

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
$\endgroup$

migrated from stackoverflow.com Apr 24 '16 at 18:55

This question came from our site for professional and enthusiast programmers.

  • 2
    $\begingroup$ 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 $\endgroup$ – user20650 Apr 18 '16 at 17:38
  • $\begingroup$ 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. $\endgroup$ – Juanchi Apr 18 '16 at 18:30
  • 1
    $\begingroup$ 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. $\endgroup$ – user20650 Apr 18 '16 at 18:48
  • $\begingroup$ 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.) $\endgroup$ – jaimedash Apr 19 '16 at 0:27
  • $\begingroup$ second choice @jaimedash ... optimum temperature for each particular strains. Edited question. Thanks. $\endgroup$ – Juanchi Apr 19 '16 at 12:50
2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Right, but as I said in my question, I have 20 strains and I just gave 3 of them for the reproducible example.. $\endgroup$ – Juanchi Apr 22 '16 at 12:37
  • $\begingroup$ I hadn't absorbed that part of the question, clearly $\endgroup$ – jaimedash Apr 22 '16 at 19:14
2
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Juanchi Apr 20 '16 at 10:56
  • $\begingroup$ 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' $\endgroup$ – Juanchi Apr 20 '16 at 11:02
  • $\begingroup$ nice. augPred is a useful command. @Juanchi Also FWIW I get no errors from augPred after running the commands here $\endgroup$ – jaimedash Apr 21 '16 at 0:26
  • $\begingroup$ 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) $\endgroup$ – jaimedash Apr 21 '16 at 0:32
  • $\begingroup$ Individual fixed model would be the suitable approach (edited answer detailing it). Thanks @jaimedash and @Walmes! $\endgroup$ – Juanchi Apr 22 '16 at 12:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.