Modelos lineales mixtos con correlacion espacial.

library(dplyr)
library(ggplot2)

library(nlme)
library(emmeans)

Ensayo evaluación de cultivares de trigo

  • Se evaluaron 15 cultivares comerciales y 35 líneas experimentales de trigo.
  • El ensayo se realizó bajo un DBCA con 3 repeticiones.
  • Se registró la posición espacial de cada parcela, y columnas (x) y filas (y).

Estadistica descriptiva

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

datos_trigo <- 
  datos_trigo |>
  mutate(Bloque = as.character(Bloque))
head(datos_trigo)
Genotipo Columna1 Bloque_vis Bloque Fila Rendimiento
G01 1 1 1 1 3.37917
G02 2 1 1 2 3.58748
G03 3 1 1 3 4.32812
G04 4 1 1 4 4.53642
G05 5 1 1 5 2.59224
G06 6 1 1 6 2.40708 
ggplot(datos_trigo, aes(Bloque, Fila)) +
  geom_raster(aes(fill = Genotipo))
ggplot(datos_trigo, aes(Bloque, Fila)) +
  geom_raster(aes(fill = Rendimiento))

ggplot(datos_trigo, aes(Genotipo, Rendimiento)) +
  stat_summary(fun.data = mean_se)

Ajuste de modelo

Sin contemplar correlación espacial

mdl_errores_indep <- lme(Rendimiento ~ Genotipo, 
list(Bloque=pdIdent(~1)),
data = datos_trigo)
datos_trigo$resid_indep <- 
  residuals(mdl_errores_indep, type = 'normalized')
datos_trigo$pred_indep <- 
  predict(mdl_errores_indep)

ggplot(datos_trigo, aes(sample = resid_indep)) +
        geom_qq(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        geom_qq_line() +
        theme_bw() +
        labs(y = "Standardised Residual", 
        x = "Theoretical")
ggplot(
        data = datos_trigo,
        mapping = aes(x = pred_indep, y = resid_indep)
        ) +
        geom_point(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        theme_bw() +
        labs(y = "Standardised Residual",
         x = "Fitted Value")
ggplot(
        data = datos_trigo,
        mapping = aes(x = Bloque, y = resid_indep)
        ) +
        geom_boxplot(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Bloque")
ggplot(
        data = datos_trigo,
        mapping = aes(x = Genotipo, y = resid_indep)
        ) +
        geom_boxplot(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Genotipo")

joint_tests(mdl_errores_indep)
model term df1 df2 F.ratio p.value
Genotipo 49 98 1.809 0.0066113 
emmeans(mdl_errores_indep, 'Genotipo') |>
  as.data.frame() |>
  arrange(emmean)
Genotipo emmean SE df lower.CL upper.CL
G36 1.465850 0.4150057 2 -0.3197755 3.251476
G30 1.604720 0.4150057 2 -0.1809055 3.390346
G33 1.712730 0.4150057 2 -0.0728955 3.498356
G45 1.743590 0.4150057 2 -0.0420355 3.529215
G46 1.759020 0.4150057 2 -0.0266055 3.544646
G28 1.774450 0.4150057 2 -0.0111755 3.560075
G38 1.789880 0.4150057 2 0.0042545 3.575505
G48 1.789880 0.4150057 2 0.0042545 3.575505
G49 1.789880 0.4150057 2 0.0042545 3.575505
G27 1.805310 0.4150057 2 0.0196845 3.590936
G23 1.851600 0.4150057 2 0.0659745 3.637226
G35 1.990470 0.4150057 2 0.2048445 3.776095
G37 2.044477 0.4150057 2 0.2588511 3.830102
G19 2.083050 0.4150057 2 0.2974245 3.868676
G29 2.083050 0.4150057 2 0.2974245 3.868676
G24 2.083050 0.4150057 2 0.2974245 3.868676
G34 2.113910 0.4150057 2 0.3282845 3.899535
G47 2.221920 0.4150057 2 0.4362945 4.007545
G25 2.252780 0.4150057 2 0.4671545 4.038405
G13 2.252780 0.4150057 2 0.4671545 4.038405
G07 2.268210 0.4150057 2 0.4825845 4.053835
G05 2.299070 0.4150057 2 0.5134445 4.084695
G21 2.329930 0.4150057 2 0.5443045 4.115556
G22 2.345360 0.4150057 2 0.5597345 4.130986
G43 2.399367 0.4150057 2 0.6137411 4.184992
G26 2.422510 0.4150057 2 0.6368845 4.208136
G31 2.445657 0.4150057 2 0.6600311 4.231282
G06 2.530520 0.4150057 2 0.7448945 4.316146
G40 2.530520 0.4150057 2 0.7448945 4.316146
G44 2.561380 0.4150057 2 0.7757545 4.347005
G15 2.561380 0.4150057 2 0.7757545 4.347005
G10 2.592240 0.4150057 2 0.8066145 4.377866
G39 2.638530 0.4150057 2 0.8529045 4.424156
G41 2.669390 0.4150057 2 0.8837645 4.455016
G50 2.669390 0.4150057 2 0.8837645 4.455016
G14 2.731110 0.4150057 2 0.9454845 4.516736
G08 2.754257 0.4150057 2 0.9686311 4.539882
G16 2.893127 0.4150057 2 1.1075011 4.678752
G42 2.900840 0.4150057 2 1.1152145 4.686465
G12 3.039710 0.4150057 2 1.2540845 4.825336
G03 3.047427 0.4150057 2 1.2618011 4.833052
G18 3.062857 0.4150057 2 1.2772311 4.848482
G02 3.093717 0.4150057 2 1.3080911 4.879342
G09 3.224870 0.4150057 2 1.4392445 5.010496
G11 3.255730 0.4150057 2 1.4701045 5.041355
G20 3.263447 0.4150057 2 1.4778211 5.049072
G32 3.294307 0.4150057 2 1.5086811 5.079932
G01 3.402317 0.4150057 2 1.6166911 5.187942
G04 3.564330 0.4150057 2 1.7787045 5.349956
G17 3.610620 0.4150057 2 1.8249945 5.396246

Contemplando correlación espacial

mdl_errores_corr_exp <- lme(Rendimiento ~ Genotipo, 
random = list(Bloque=pdIdent(~1)),
correlation = corExp(form=~Bloque_vis + Fila,
nugget=FALSE,
metric='euclidean'),
data = datos_trigo)
mdl_errores_corr_gau <- lme(Rendimiento ~ Genotipo, 
random = list(Bloque=pdIdent(~1)),
correlation = corGaus(form=~Bloque_vis + Fila,
nugget=FALSE,
metric='euclidean'),
data = datos_trigo)

mdl_errores_corr_lin <- lme(Rendimiento ~ Genotipo, 
random = list(Bloque=pdIdent(~1)),
correlation = corLin(form=~Bloque_vis + Fila,
nugget=FALSE,
metric='euclidean'),
data = datos_trigo)

mdl_errores_corr_spher <- lme(Rendimiento ~ Genotipo, 
random = list(Bloque=pdIdent(~1)),
correlation = corSpher(form=~Bloque_vis + Fila,
nugget=FALSE,
metric='euclidean'),
data = datos_trigo)

mdl_errores_corr_rationalquad<- lme(Rendimiento ~ Genotipo, 
random = list(Bloque=pdIdent(~1)),
correlation = corRatio(form=~Bloque_vis + Fila,
nugget=FALSE,
metric='euclidean'),
data = datos_trigo)

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 (Linear) 305.3934 443.4674
Correlacion Espacial (Exp) 306.1659 444.2400
Correlacion Espacial (Spheric) 311.3027 449.3767
Correlacion Espacial (Rational Quadratic) 316.7885 454.8625
Correlacion Espacial (Gaussian) 334.8265 472.9005
Errores Independientes 376.6707 512.1395

Modelo final

datos_trigo$resid_corr_lin <- 
  residuals(mdl_errores_corr_lin, type = 'normalized')
datos_trigo$pred_corr_lin <- 
  predict(mdl_errores_corr_lin)

ggplot(datos_trigo, aes(sample = resid_corr_lin)) +
        geom_qq(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        geom_qq_line() +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Theoretical")
ggplot(
        data = datos_trigo,
        mapping = aes(x = pred_corr_lin, y = resid_corr_lin)
        ) +
        geom_point(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Fitted Value")
ggplot(
        data = datos_trigo,
        mapping = aes(x = Bloque, y = resid_corr_lin)
        ) +
        geom_boxplot(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Bloque")
ggplot(
        data = datos_trigo,
        mapping = aes(x = Genotipo, y = resid_corr_lin)
        ) +
        geom_boxplot(
            colour = "black",
            fill = "aquamarine3",
            size = 2,
            shape = 21
        ) +
        theme_bw() +
        labs(y = "Standardised Residual", x = "Genotipo")

joint_tests(mdl_errores_corr_lin)
model term df1 df2 F.ratio p.value
Genotipo 49 98 6.667 0 
emmeans(mdl_errores_corr_lin, 'Genotipo') |>
  as.data.frame() |>
  arrange(emmean)
Genotipo emmean SE df lower.CL upper.CL
G28 1.476768 0.3023253 2 0.1759668 2.777568
G36 1.538156 0.2970479 2 0.2600621 2.816250
G21 1.761386 0.2940460 2 0.4962077 3.026563
G43 1.809996 0.2994555 2 0.5215430 3.098449
G46 1.840904 0.2967700 2 0.5640060 3.117802
G45 1.887223 0.3024035 2 0.5860863 3.188360
G07 1.909751 0.2942985 2 0.6434873 3.176016
G29 1.978285 0.2945636 2 0.7108797 3.245689
G30 1.998038 0.2949078 2 0.7291527 3.266924
G49 2.017462 0.2937601 2 0.7535146 3.281410
G06 2.068855 0.3003842 2 0.7764063 3.361304
G23 2.081663 0.2999518 2 0.7910741 3.372251
G38 2.136192 0.2952369 2 0.8658903 3.406494
G33 2.141480 0.2964796 2 0.8658307 3.417129
G15 2.259164 0.2972114 2 0.9803666 3.537961
G44 2.302427 0.2991998 2 1.0150747 3.589780
G31 2.326669 0.2939662 2 1.0618348 3.591504
G27 2.347612 0.3013353 2 1.0510702 3.644153
G37 2.357334 0.2973078 2 1.0781211 3.636546
G34 2.378812 0.2969033 2 1.1013406 3.656284
G48 2.381430 0.3006979 2 1.0876319 3.675229
G39 2.410992 0.3006901 2 1.1172271 3.704757
G19 2.411993 0.2896015 2 1.1659380 3.658047
G05 2.417046 0.2976442 2 1.1363868 3.697706
G25 2.419023 0.2950482 2 1.1495331 3.688513
G24 2.430671 0.2993659 2 1.1426037 3.718739
G35 2.471400 0.2925652 2 1.2125933 3.730206
G41 2.492912 0.3012054 2 1.1969293 3.788894
G47 2.522935 0.2955686 2 1.2512058 3.794664
G18 2.534374 0.3004206 2 1.2417685 3.826980
G14 2.585068 0.3022517 2 1.2845836 3.885552
G26 2.589415 0.2948478 2 1.3207879 3.858043
G08 2.601926 0.2995628 2 1.3130110 3.890840
G13 2.609047 0.2948342 2 1.3404775 3.877616
G40 2.658114 0.3044837 2 1.3480260 3.968202
G22 2.672658 0.3003799 2 1.3802276 3.965089
G16 2.797675 0.3067174 2 1.4779770 4.117374
G42 2.850870 0.2964127 2 1.5755091 4.126231
G03 2.874526 0.2945488 2 1.6071851 4.141868
G02 2.892058 0.2959940 2 1.6184990 4.165617
G10 2.899056 0.2959095 2 1.6258600 4.172251
G50 3.019208 0.2984731 2 1.7349821 4.303435
G20 3.041012 0.2950525 2 1.7715036 4.310521
G11 3.069149 0.3044282 2 1.7593005 4.378998
G32 3.069396 0.2936543 2 1.8059038 4.332889
G01 3.107186 0.3103365 2 1.7719153 4.442456
G09 3.107770 0.2973122 2 1.8285389 4.387001
G17 3.207479 0.2965052 2 1.9317205 4.483238
G04 3.363593 0.2989042 2 2.0775123 4.649674
G12 3.403372 0.2990477 2 2.1166740 4.690071