Modelos lineales mixtos con correlacion espacial.

library(dplyr)
library(ggplot2)

library(nlme)
library(emmeans)

Ensayo testigos a la par

Se evalúan 16 híbridos y un testigo común intercalado
El testigo permite modelar los efectos sistemáticos de la calidad del terreno donde se ubican las UE.

Estadistica descriptiva

datos_testigosap <- read.table('Datos/Dia 1/Testigos_apareados.txt',
sep = '\t',
header = TRUE)

head(datos_testigosap)

Posicion	Hibrido	Rendimiento	Tes	Hib	Dif	PRED_Dif
1	Testigo	896	896	886	-10	75.03804
2	H1	886	NA	NA	NA	64.55978
3	H2	1037	NA	NA	NA	54.08152
4	Testigo	1101	1101	1037	-64	43.60326
5	Testigo	1055	1055	1222	167	33.12500
6	H3	1222	NA	NA	NA	22.64674

ggplot(datos_testigosap, aes(Posicion, Rendimiento, color = Hibrido == 'Testigo')) +
  geom_point()

Ajuste de modelo

Errores independientes

mdl_errores_indep <- aov(Rendimiento ~ Hibrido, 
data = datos_testigosap)

Contemplando correlación espacial

mdl_errores_corr_exp <- gls(Rendimiento ~ Hibrido, 
correlation = corExp(form=~Posicion,
nugget=FALSE,
metric='euclidean'),
data = datos_testigosap)

mdl_errores_corr_gau <- gls(Rendimiento ~ Hibrido, 
correlation = corGaus(form=~Posicion,
nugget=FALSE,
metric='euclidean'),
data = datos_testigosap)

mdl_errores_corr_lin <- gls(Rendimiento ~ Hibrido, 
correlation = corLin(form=~Posicion,
nugget=FALSE,
metric='euclidean'),
data = datos_testigosap)

mdl_errores_corr_spher <- gls(Rendimiento ~ Hibrido, 
correlation = corSpher(form=~Posicion,
nugget=FALSE,
metric='euclidean'),
data = datos_testigosap)

mdl_errores_corr_rationalquad<- gls(Rendimiento ~ Hibrido, 
correlation = corRatio(form=~Posicion,
nugget=FALSE,
metric='euclidean'),
data = datos_testigosap)

Comparacion de modelos

data.frame(
  Modelo = c('Errores Independientes',
             'Correlacion Espacial (Exp)',
             'Correlacion Espacial (Gaussian)',
             'Correlacion Espacial (Linear)',
             'Correlacion Espacial (Spheric)',
             'Correlacion Espacial (Rational Quadratic)'),
  AIC = c(AIC(mdl_errores_indep), 
          AIC(mdl_errores_corr_exp),
          AIC(mdl_errores_corr_gau),
          AIC(mdl_errores_corr_lin),
          AIC(mdl_errores_corr_spher),
          AIC(mdl_errores_corr_rationalquad)
          ),
  BIC = c(BIC(mdl_errores_indep), 
          BIC(mdl_errores_corr_exp),
          BIC(mdl_errores_corr_gau),
          BIC(mdl_errores_corr_lin),
          BIC(mdl_errores_corr_spher),
          BIC(mdl_errores_corr_rationalquad))
) |> 
  arrange(AIC, BIC)

Modelo	AIC	BIC
Correlacion Espacial (Exp)	218.6227	232.0756
Correlacion Espacial (Rational Quadratic)	218.8059	232.2589
Correlacion Espacial (Linear)	219.1268	232.5798
Correlacion Espacial (Gaussian)	219.1685	232.6214
Correlacion Espacial (Spheric)	219.2063	232.6593
Errores Independientes	398.1584	424.5416

Modelo final

datos_testigosap$resid_corr_exp <- 
  residuals(mdl_errores_corr_exp, type = 'normalized')
datos_testigosap$pred_corr_exp <-
  predict(mdl_errores_corr_exp)

ggplot(datos_testigosap, aes(sample = resid_corr_exp)) +
        geom_qq(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        geom_qq_line() +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Theoretical")
ggplot(
        data = datos_testigosap,
        mapping = aes(x = pred_corr_exp, y = resid_corr_exp)
        ) +
        geom_point(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Fitted Value")
ggplot(
        data = datos_testigosap,
        mapping = aes(x = Hibrido, y = resid_corr_exp)
        ) +
        geom_boxplot(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Hibrido")

joint_tests(mdl_errores_corr_exp)

model term	df1	df2	F.ratio	p.value
Hibrido	16	6.9	5.273	0.0171447

emmeans(mdl_errores_corr_exp, 'Hibrido') |>
  as.data.frame() |>
  arrange(emmean)

Hibrido	emmean	SE	df	lower.CL	upper.CL
H15	653.3557	85.33059	13.968702	470.3013	836.4101
H13	727.5479	85.33059	13.968702	544.4935	910.6023
H16	902.8669	85.33059	13.968702	719.8125	1085.9213
H8	933.4034	85.33059	13.968702	750.3490	1116.4577
H7	936.2055	85.33059	13.968702	753.1511	1119.2599
H6	966.3110	85.33059	13.968702	783.2566	1149.3654
H12	979.7981	85.33059	13.968702	796.7437	1162.8525
H1	1005.4633	85.33059	13.968702	822.4089	1188.5177
H14	1070.0660	85.33059	13.968702	887.0116	1253.1204
H9	1078.2818	85.33059	13.968702	895.2274	1261.3362
H11	1078.4279	85.33059	13.968702	895.3736	1261.4823
H2	1091.0732	85.33059	13.968702	908.0188	1274.1276
Testigo	1096.9845	45.08659	1.493844	824.1424	1369.8265
H5	1128.6462	85.33059	13.968702	945.5918	1311.7006
H10	1145.6426	85.33059	13.968702	962.5882	1328.6970
H4	1244.1896	85.33059	13.968702	1061.1352	1427.2440
H3	1248.3107	85.33059	13.968702	1065.2563	1431.3651