No regresión lineal modelo mixto

Question

Tengo un datastet con varias cepas de hongos (20) y estoy tratando de averiguar cual es la temperatura óptima de crecimiento de cada cepa. Para eso, les pongo el crecimiento en las placas a diferentes temperaturas (22,25,28 y 31ºC).
A continuación es de un reducido conjunto de datos con 3 sitios:

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)

Yo se ajustan a un modelo con un único nivel de factor de

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

Hay una manera de (modelos mixtos?) para adaptarse a un modelo único de cuentas para todas las cepas de hongos?

Como ya he dicho, tengo particular interés en la thx parámetro (temperatura Óptima) y su intervalo de confianza. Podría alguien ayudarme con eso?

Tras la respuesta de comentario tengo los parámetros de thx, y estoy tratando de conseguir esta tabla de salida:

#       thx     lwr   upr
# Iti   25.13   -     -
# Lon   24.44   -     -
# Sinop 26.33   -     -

Este es un Walmes solución (bandeja de entrada) que coincide con @jaimedash uno

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

Answer 1

Porque desea predecir las temperaturas óptimas de cada cepa, el tratamiento de las cepas como fijo tiene sentido. Sin embargo, una interacción entre los tres niveles del factor y el modelo para la óptima hace que para un dolor de cabeza. Creo que no hay suficientes datos para ajustarlo.

Un modelo de efectos aleatorios podría ayudar. Con sus datos completos, de 20 de grupos a los que usted puede ser capaz de adaptarse fácilmente a una gran efectos aleatorios estructura como en la otra respuesta (o incluso mayor). Pero, en el caso de que usted no recibe mucha ayuda a la predicción específica óptima de cada cepa. (A pesar de que había estimado la media de los parámetros del modelo óptimo de forma más precisa debido a la puesta en común de datos entre los tres grupos.

El enfoque más sencillo es ajustar un modelo independiente para cada cepa y predecir, a partir de eso. Aquí está el código adoptado a partir de la tuya (y también el uso de intervals como en @Walmes respuesta) que hace eso y combina las estimaciones para obtener la tabla de estimaciones e intervalos de tiempo:

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

tabla:

      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

usted puede explorar el ajuste de cada modelo de uso de la plot(reslist[[1]]$model) o summary(reslist[[1]])

editar oh, como @Walmes señalado (y la pregunta es actualmente editado para mostrar) usted puede hacer esto de una manera más limpia usando nlme:nlsList