Verificación de supuestos

Se utilizó la libraría DHARMa para el análisis de residuos. Esta librería utiliza un método basado en simulación para reescalar los residuos entre 0 y 1 a partir de los cuantiles.

Con estos residuos se evidenciaron problemas en los supuestos (e.g., distribución, dispersión).

# Parámetro de Dispersion 

e1 <- resid(m1.cu, type = "pearson")
n  <- nrow(cobre)
k  <- length(fixef(m1.cu)) + 1 #parametros estimados; se suma 1varianza(Num replicas)
dispersion <- sum(e1^2) / (n - k)

Por otro lado, la estimación del parámetro de dispersión arrojó un valor de 0.236, por lo que se considera que el modelo presenta subdispersión.

Verificación de supuestos

Se volvió a utilizar la libraría DHARMa para el análisis de residuos.

A partir del análisis gráfico de los residuos de este modelo también se evidenciaron problemas en los supuestos.

#  Parámetro de Dispersion 

e1 <- resid(m1.cu_basal, type = "pearson")
n  <- nrow(cobre_medrep)
k  <- length(fixef(m1.cu_basal)) + 1 #parametros estimados; se suma 1varianza(Num replicas)
dispersion <- sum(e1^2) / (n - k)

Por otro lado, la estimación del parámetro de dispersión arrojó un valor de 0.419, por lo que mejora respecta del modelo anterior, pero no demasiado.

La diferencia en el AIC entre el primer modelo (446.48) y el segundo (356.09) es de 90.39, por lo que el modelo con el nivel basal como predictor ajusta mejor.

Verificación de supuestos

# Con DHARMA

simulationOutput <- simulateResiduals(m2.cu)
plot(simulationOutput)

hist(simulationOutput)

e1 <- resid(m2.cu, type = "pearson")
n  <- nrow(cobre)
k  <- length(fixef(m2.cu)) + 1 #parametros estimados; se suma 1varianza(Num replicas)
dispersion <- sum(e1^2) / (n - k)

# Con DHARMA

simulationOutput <- simulateResiduals(m2.cu_basal)
plot(simulationOutput)

hist(simulationOutput)

e1 <- resid(m2.cu_basal, type = "pearson")
n  <- nrow(cobre_medrep)
k  <- length(fixef(m2.cu_basal)) + 1 #parametros estimados; se suma 1varianza(Num replicas)
dispersion <- sum(e1^2) / (n - k)

Selección de modelo por LOOCV

Se calculó el error cuadrático medio mediante la técnica Leave-one-out cross-validation para los dos últimos modelos y así seleccionar el modelo que prediga de forma más robusta la cantidad de hojas por réplica a partir del tratamiento y la semana. Para ello, se ajustó el modelo con todas las observaciones excepto una y se comparó la cantidad de hojas en ese tratamiento en esa semana con la predicción realizada a partir del modelo. Este proceso se repitió para todos tratamientos y todas las semanas. Luego se promedió este error cuadrático. El modelo que mejor ajuste será aquel que tenga un menor error cuadrático medio.

# cinco semanas
ec <- c()

for(i in levels(droplevels(cobre$Replica))){
  for(j in levels(droplevels(cobre$Semana))){
    
    cobre.loocv <- glmmTMB(Nhojas ~ Tratamiento*Semana + (1|Replica), data = cobre[cobre$Replica!=i|cobre$Semana!=j,], family = compois, REML = TRUE)
    
    ec <- c(ec,(predict(cobre.loocv, newdata = cobre[cobre$Replica==i&cobre$Semana==j,])-cobre$Nhojas[cobre$Replica==i&cobre$Semana==j])^2)

  }
}

ecm <- mean(ec)




# sacando basal
cobre_medrep$Semana <- as.factor(cobre_medrep$Semana)

ec_basal <- c()

for(i in levels(droplevels(cobre_medrep$Replica))){
  for(j in levels(droplevels(cobre_medrep$Semana))){
    
    cobre.loocv <- glmmTMB(Nhojas ~ Tratamiento*Semana+Basal + (1|Replica), data = cobre_medrep[cobre_medrep$Replica!=i|cobre_medrep$Semana!=j,], family = compois, REML = TRUE)
    
    ec_basal <- c(ec_basal,(predict(cobre.loocv, newdata = cobre_medrep[cobre_medrep$Replica==i&cobre_medrep$Semana==j,])-cobre_medrep$Nhojas[cobre_medrep$Replica==i&cobre_medrep$Semana==j])^2)

  }
}


ecm_basal <- mean(ec_basal)

El error cuadrático medio del modelo sin controlar por el nivel basal de hojas por réplica (ECM = 552.41) fue menor que aquel obtenido controlando por dicho nivel (ECM = 645.3).

Conclusión: Dados estos análisis, se decide avanzar con el modelo de Conway y Maxwell sin controlar por el nivel basal.

Gráficos

# Medias con IC
gfxcomp.cu <-as.data.frame(comp2$emmeans)


ggplot(gfxcomp.cu, aes(Tratamiento, y=response)) + 
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL,colour = factor(gfxcomp.cu$Semana)), width=.1) +
  geom_point(aes(colour = factor(gfxcomp.cu$Semana)), size = 4)+ ylab("Numero de hojas predichas por replica") + 
  xlab("Tratamiento") +
  theme_grey(base_size = 10) +  
  ggtitle("Número de hojas predichas según semana en cada tratamiento",
          "Media (EE)") + guides(color=guide_legend("Numero de hojas\nen c/semana"))+
  theme(axis.text.x  = element_text(vjust=0.5, size=10,angle=0, ))+ 
  scale_color_manual(values=c("#264653", "#2a9d8f", "#e9c46a", "#f4a261", "#e76f51"))

ggplot(gfxcomp.cu, aes(Semana, y=response)) + 
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL,colour = factor(gfxcomp.cu$Tratamiento)), width=.1) +
  geom_point(aes(colour = factor(gfxcomp.cu$Tratamiento)), size = 4)+ ylab("Numero de hojas predichas por replica") + 
  xlab("Semana") +
  theme_grey(base_size = 16) +  
  ggtitle("Número de hojas predichas según tratamientos en cada semana",
          "Media (EE)") + guides(color=guide_legend("Numero de hojas\nsegún tratamiento"))+
  theme(plot.title = element_text(size = 12),
        plot.subtitle = element_text(size = 10),
        axis.title = element_text(size=12),
        legend.title = element_text(size=12),
    axis.text.x  = element_text(vjust=0.5, size=10))+ 
  scale_color_manual(values=c("#264653", "#e9c46a", "#e76f51"))

Verificación de supuestos

También para los modelos de Zinc se utilizó la libraría DHARMa para el análisis de residuos.

Con estos residuos también se evidenciaron problemas en los supuestos.

# dispersion
e1 <- resid(m1.zn, type = "pearson")
n  <- nrow(zinc)
k  <- length(fixef(m1.zn)) + 1 #parametros estimados; se suma 1varianza(Num replicas)
dispersion <- sum(e1^2) / (n - k)

Por otro lado, la estimación del parámetro de dispersión arrojó un valor de 0.347, por lo que se considera que el modelo presenta subdispersión.

Verificación de supuestos

Nuevamente, los residuos a partir de simulaciones presentaron problemas.

# Dispersion 

e1 <- resid(m1.zn_basal, type = "pearson")
n  <- nrow(zinc_medrep)
k  <- length(fixef(m1.zn_basal)) + 1 #parametros estimados; se suma 1varianza(Num replicas)
dispersion <- sum(e1^2) / (n - k)

Por otro lado, la estimación del parámetro de dispersión arrojó un valor de 0.33, por lo que empeora respecta del modelo anterior, pero no demasiado.

La diferencia en el AIC entre el primer modelo (447) y el segundo (369.74) es de 77.25, por lo que el modelo con el nivel basal como predictor ajusta mejor.

Verificación de supuestos

# Con DHARMA

simulationOutput <- simulateResiduals(m2.zn)
plot(simulationOutput)

hist(simulationOutput)

e1 <- resid(m2.zn, type = "pearson")
n  <- nrow(zinc)
k  <- length(fixef(m2.zn)) + 1 #parametros estimados; se suma 1varianza(Num replicas)
dispersion <- sum(e1^2) / (n - k)

# Con DHARMA

simulationOutput <- simulateResiduals(m2.zn_basal)
plot(simulationOutput)

hist(simulationOutput)

e1 <- resid(m2.zn_basal, type = "pearson")
n  <- nrow(zinc_medrep)
k  <- length(fixef(m2.zn_basal)) + 1 #parametros estimados; se suma 1varianza(Num replicas)
dispersion <- sum(e1^2) / (n - k)

Selección de modelo por LOOCV

Se calculó el error cuadrático medio mediante la técnica Leave-one-out cross-validation para los dos últimos modelos y así seleccionar el modelo que prediga de forma más robusta la cantidad de hojas por réplica a partir del tratamiento y la semana.

# cinco semanas
ec <- c()

for(i in levels(droplevels(zinc$Replica))){
  for(j in levels(droplevels(zinc$Semana))){
    
    zinc.loocv <- glmmTMB(Nhojas ~ Tratamiento*Semana + (1|Replica), data = zinc[zinc$Replica!=i|zinc$Semana!=j,], family = compois, REML = TRUE)
    
    ec <- c(ec,(predict(zinc.loocv, newdata = zinc[zinc$Replica==i&zinc$Semana==j,])-zinc$Nhojas[zinc$Replica==i&zinc$Semana==j])^2)

  }
}

ecm <- mean(ec)

# sacando basal
zinc_medrep$Semana <- as.factor(zinc_medrep$Semana)

ec_basal <- c()

for(i in levels(droplevels(zinc_medrep$Replica))){
  for(j in levels(droplevels(zinc_medrep$Semana))){
    
    zinc.loocv <- glmmTMB(Nhojas ~ Tratamiento*Semana+Basal + (1|Replica), data = zinc_medrep[zinc_medrep$Replica!=i|zinc_medrep$Semana!=j,], family = compois, REML = TRUE)
    
    ec_basal <- c(ec_basal,(predict(zinc.loocv, newdata = zinc_medrep[zinc_medrep$Replica==i&zinc_medrep$Semana==j,])-zinc_medrep$Nhojas[zinc_medrep$Replica==i&zinc_medrep$Semana==j])^2)

  }
}


ecm_basal <- mean(ec_basal)

El error cuadrático medio del modelo sin controlar por el nivel basal de hojas por réplica (ECM = 691.73) fue menor que aquel obtenido controlando por dicho nivel (ECM = 820.73).

Conclusión: Dados estos análisis, se decide avanzar con el modelo de Conway y Maxwell sin controlar por el nivel basal.

Gráficos

# Medias con IC
gfxcomp.zn <-as.data.frame(comp2$emmeans)


ggplot(gfxcomp.zn, aes(factor(Tratamiento, levels = c("Control","Zn500","Zn1700")), y=response)) + 
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL,colour = factor(gfxcomp.zn$Semana)), width=.1) +
  geom_point(aes(colour = factor(gfxcomp.zn$Semana)), size = 4)+ ylab("Numero de hojas predichas por replica") + 
  xlab("Tratamiento") +
  theme_grey(base_size = 10) +  ggtitle("Número de hojas predichas según semana en cada tratamiento","Media (EE)") + guides(color=guide_legend("Numero de hojas\nen cada semana"))+
  theme(axis.text.x  = element_text(vjust=0.5, size=12,angle=0))+ 
  scale_color_manual(values=c("#264653", "#2a9d8f", "#e9c46a", "#f4a261", "#e76f51"))

ggplot(gfxcomp.zn, aes(Semana, y=response)) + 
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL,colour = factor(gfxcomp.zn$Tratamiento, levels = c("Control","Zn500","Zn1700"))), width=.1) +
  geom_point(aes(colour = factor(gfxcomp.zn$Tratamiento, levels = c("Control","Zn500","Zn1700"))), size = 4)+ ylab("Numero de hojas predichas por replica") + 
  xlab("Semana") +
  theme_grey(base_size = 16) +  ggtitle("Número de hojas predichas según tratamientos en cada semana", "Media (EE)") + guides(color=guide_legend("Numero de hojas\nsegún tratamiento"))+
  theme(plot.title = element_text(size = 12),
        plot.subtitle = element_text(size = 10),
        axis.title = element_text(size=12),
        legend.title = element_text(size=12),
        axis.text.x  = element_text(vjust=0.5, size=10))+ 
  scale_color_manual(values=c("#264653", "#e9c46a", "#e76f51"))