setwd("G:/Mi unidad/Agrosavia/FeCa/Fenoma/Análisis/AguaS")
datos<-read.table("humecta.csv", header=T, sep=';')
##Librerias
library(lme4)
## Cargando paquete requerido: Matrix
library(lmerTest) # p-values
## Warning: package 'lmerTest' was built under R version 4.4.3
##
## Adjuntando el paquete: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(emmeans) # post hoc
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(dplyr)
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
##
## Adjuntando el paquete: 'tidyr'
## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack
library(agricolae)
## Warning: package 'agricolae' was built under R version 4.4.3
library(metan)
## Warning: package 'metan' was built under R version 4.4.3
## |=========================================================|
## | Multi-Environment Trial Analysis (metan) v1.19.0 |
## | Author: Tiago Olivoto |
## | Type 'citation('metan')' to know how to cite metan |
## | Type 'vignette('metan_start')' for a short tutorial |
## | Visit 'https://bit.ly/metanpkg' for a complete tutorial |
## |=========================================================|
##
## Adjuntando el paquete: 'metan'
## The following object is masked from 'package:tidyr':
##
## replace_na
## The following object is masked from 'package:dplyr':
##
## recode_factor
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
##Convertir a factor
datos$gen <- factor(datos$gen)
datos$municipio <- factor(datos$municipio)
datos$mun <- factor(datos$mun)
datos$reg <- factor(datos$reg)
#Modelo 0
modelo <- lm (log(teta) ~ gen + mun,
data = datos)
anova(modelo)
## Analysis of Variance Table
##
## Response: log(teta)
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 7 2.48313 0.35473 1623.079 < 2.2e-16 ***
## mun 9 0.04800 0.00533 24.403 < 2.2e-16 ***
## Residuals 60 0.01311 0.00022
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Modelo bruto
modelo <- lm (teta ~ gen + mun,
data = datos)
anova(modelo)
## Analysis of Variance Table
##
## Response: teta
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 7 12219.6 1745.66 2605.387 < 2.2e-16 ***
## mun 9 270.3 30.04 44.828 < 2.2e-16 ***
## Residuals 60 40.2 0.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##Contrastes a posteriori
#Genotipos
g<-emmeans(modelo, pairwise ~ gen)
g
## $emmeans
## gen emmean SE df lower.CL upper.CL
## CNCH12 86.2 0.259 60 85.7 86.7
## CNCH13 88.5 0.259 60 88.0 89.0
## FBO1 72.6 0.259 60 72.1 73.1
## FCHI8 94.7 0.294 60 94.2 95.3
## FEAR5 86.3 0.259 60 85.8 86.8
## FGI4 51.9 0.259 60 51.4 52.4
## FMA7 77.0 0.259 60 76.5 77.5
## FSV1 86.9 0.275 60 86.4 87.5
##
## Results are averaged over the levels of: mun
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## CNCH12 - CNCH13 -2.280 0.366 60 -6.228 <.0001
## CNCH12 - FBO1 13.580 0.366 60 37.097 <.0001
## CNCH12 - FCHI8 -8.529 0.392 60 -21.783 <.0001
## CNCH12 - FEAR5 -0.110 0.366 60 -0.300 1.0000
## CNCH12 - FGI4 34.310 0.366 60 93.726 <.0001
## CNCH12 - FMA7 9.220 0.366 60 25.187 <.0001
## CNCH12 - FSV1 -0.713 0.378 60 -1.888 0.5640
## CNCH13 - FBO1 15.860 0.366 60 43.326 <.0001
## CNCH13 - FCHI8 -6.249 0.392 60 -15.960 <.0001
## CNCH13 - FEAR5 2.170 0.366 60 5.928 <.0001
## CNCH13 - FGI4 36.590 0.366 60 99.955 <.0001
## CNCH13 - FMA7 11.500 0.366 60 31.415 <.0001
## CNCH13 - FSV1 1.567 0.378 60 4.147 0.0026
## FBO1 - FCHI8 -22.109 0.392 60 -56.464 <.0001
## FBO1 - FEAR5 -13.690 0.366 60 -37.398 <.0001
## FBO1 - FGI4 20.730 0.366 60 56.629 <.0001
## FBO1 - FMA7 -4.360 0.366 60 -11.910 <.0001
## FBO1 - FSV1 -14.293 0.378 60 -37.839 <.0001
## FCHI8 - FEAR5 8.419 0.392 60 21.502 <.0001
## FCHI8 - FGI4 42.839 0.392 60 109.405 <.0001
## FCHI8 - FMA7 17.749 0.392 60 45.329 <.0001
## FCHI8 - FSV1 7.816 0.399 60 19.569 <.0001
## FEAR5 - FGI4 34.420 0.366 60 94.027 <.0001
## FEAR5 - FMA7 9.330 0.366 60 25.487 <.0001
## FEAR5 - FSV1 -0.603 0.378 60 -1.597 0.7500
## FGI4 - FMA7 -25.090 0.366 60 -68.540 <.0001
## FGI4 - FSV1 -35.023 0.378 60 -92.717 <.0001
## FMA7 - FSV1 -9.933 0.378 60 -26.297 <.0001
##
## Results are averaged over the levels of: mun
## P value adjustment: tukey method for comparing a family of 8 estimates
pwpp(g, type = "response")
## Warning: `aes_()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`
## ℹ The deprecated feature was likely used in the emmeans package.
## Please report the issue at <https://github.com/rvlenth/emmeans/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#Municipios
m<-emmeans(modelo, pairwise ~ mun)
m
## $emmeans
## mun emmean SE df lower.CL upper.CL
## Chi 77.2 0.289 60 76.7 77.8
## Gig 81.3 0.289 60 80.7 81.9
## HtC 83.4 0.289 60 82.9 84.0
## Jam 79.8 0.339 60 79.2 80.5
## PtR 81.0 0.312 60 80.4 81.7
## RiN 79.9 0.289 60 79.3 80.5
## SnV 78.1 0.289 60 77.6 78.7
## Tam 83.0 0.289 60 82.4 83.6
## ViG 81.4 0.289 60 80.9 82.0
## Yac 80.0 0.289 60 79.4 80.5
##
## Results are averaged over the levels of: gen
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Chi - Gig -4.0375 0.409 60 -9.865 <.0001
## Chi - HtC -6.2000 0.409 60 -15.149 <.0001
## Chi - Jam -2.5978 0.446 60 -5.828 <.0001
## Chi - PtR -3.7931 0.425 60 -8.916 <.0001
## Chi - RiN -2.6750 0.409 60 -6.536 <.0001
## Chi - SnV -0.9000 0.409 60 -2.199 0.4688
## Chi - Tam -5.7625 0.409 60 -14.080 <.0001
## Chi - ViG -4.2000 0.409 60 -10.262 <.0001
## Chi - Yac -2.7125 0.409 60 -6.628 <.0001
## Gig - HtC -2.1625 0.409 60 -5.284 0.0001
## Gig - Jam 1.4397 0.446 60 3.230 0.0578
## Gig - PtR 0.2444 0.425 60 0.575 0.9999
## Gig - RiN 1.3625 0.409 60 3.329 0.0446
## Gig - SnV 3.1375 0.409 60 7.666 <.0001
## Gig - Tam -1.7250 0.409 60 -4.215 0.0032
## Gig - ViG -0.1625 0.409 60 -0.397 1.0000
## Gig - Yac 1.3250 0.409 60 3.237 0.0567
## HtC - Jam 3.6022 0.446 60 8.081 <.0001
## HtC - PtR 2.4069 0.425 60 5.658 <.0001
## HtC - RiN 3.5250 0.409 60 8.613 <.0001
## HtC - SnV 5.3000 0.409 60 12.950 <.0001
## HtC - Tam 0.4375 0.409 60 1.069 0.9859
## HtC - ViG 2.0000 0.409 60 4.887 0.0003
## HtC - Yac 3.4875 0.409 60 8.521 <.0001
## Jam - PtR -1.1952 0.457 60 -2.613 0.2343
## Jam - RiN -0.0772 0.446 60 -0.173 1.0000
## Jam - SnV 1.6978 0.446 60 3.809 0.0114
## Jam - Tam -3.1647 0.446 60 -7.100 <.0001
## Jam - ViG -1.6022 0.446 60 -3.594 0.0214
## Jam - Yac -0.1147 0.446 60 -0.257 1.0000
## PtR - RiN 1.1181 0.425 60 2.628 0.2276
## PtR - SnV 2.8931 0.425 60 6.801 <.0001
## PtR - Tam -1.9694 0.425 60 -4.630 0.0008
## PtR - ViG -0.4069 0.425 60 -0.957 0.9936
## PtR - Yac 1.0806 0.425 60 2.540 0.2691
## RiN - SnV 1.7750 0.409 60 4.337 0.0021
## RiN - Tam -3.0875 0.409 60 -7.544 <.0001
## RiN - ViG -1.5250 0.409 60 -3.726 0.0146
## RiN - Yac -0.0375 0.409 60 -0.092 1.0000
## SnV - Tam -4.8625 0.409 60 -11.881 <.0001
## SnV - ViG -3.3000 0.409 60 -8.063 <.0001
## SnV - Yac -1.8125 0.409 60 -4.429 0.0016
## Tam - ViG 1.5625 0.409 60 3.818 0.0111
## Tam - Yac 3.0500 0.409 60 7.452 <.0001
## ViG - Yac 1.4875 0.409 60 3.634 0.0190
##
## Results are averaged over the levels of: gen
## P value adjustment: tukey method for comparing a family of 10 estimates
pwpp(m, type = "response")
# Modelo 1
modelo <- lmer(log(teta) ~ gen +
(1|mun),
data = datos)
anova(modelo)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## gen 2.4675 0.3525 7 60.066 1613.3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ranova(modelo)
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log(teta) ~ gen + (1 | mun)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 10 169.63 -319.27
## (1 | mun) 9 135.55 -253.10 68.163 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Modelo 2
modelo_blup <- lmer(teta ~ 1 +
(1|gen) +
(1|mun),
data = datos)
ranova(modelo_blup)
## boundary (singular) fit: see help('isSingular')
## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## teta ~ (1 | gen) + (1 | mun)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 -139.23 286.46
## (1 | gen) 3 -304.01 614.01 329.55 1 < 2.2e-16 ***
## (1 | mun) 3 -188.02 382.03 97.57 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Extraer Blups
blups <- ranef(modelo_blup)
#Blups Gen
blups$gen
## (Intercept)
## CNCH12 5.682404
## CNCH13 7.961560
## FBO1 -7.892568
## FCHI8 14.207994
## FEAR5 5.792363
## FGI4 -28.614893
## FMA7 -3.534182
## FSV1 6.397323
#Valor predicho
fixef(modelo_blup)[1] + blups$gen
## (Intercept)
## CNCH12 86.20790
## CNCH13 88.48705
## FBO1 72.63292
## FCHI8 94.73348
## FEAR5 86.31785
## FGI4 51.91060
## FMA7 76.99131
## FSV1 86.92281
#Blups Parcela
blups$mun
## (Intercept)
## Chi -3.2150660
## Gig 0.7328851
## HtC 2.8474224
## Jam -0.6709794
## PtR 0.4914815
## RiN -0.5993956
## SnV -2.3350274
## Tam 2.4196258
## ViG 0.8917810
## Yac -0.5627273
fixef(modelo_blup)[1] + blups$mun
## (Intercept)
## Chi 77.31043
## Gig 81.25838
## HtC 83.37291
## Jam 79.85451
## PtR 81.01697
## RiN 79.92610
## SnV 78.19046
## Tam 82.94512
## ViG 81.41727
## Yac 79.96276
#Tabla blup_gen
blup_gen <- ranef(modelo_blup)$gen %>%
tibble::rownames_to_column("gen") %>%
rename(BLUP = `(Intercept)`)
blup_gen
## gen BLUP
## 1 CNCH12 5.682404
## 2 CNCH13 7.961560
## 3 FBO1 -7.892568
## 4 FCHI8 14.207994
## 5 FEAR5 5.792363
## 6 FGI4 -28.614893
## 7 FMA7 -3.534182
## 8 FSV1 6.397323
#Tabla blup_mun
blup_mun <- ranef(modelo_blup)$mun %>%
tibble::rownames_to_column("mun") %>%
rename(BLUP = `(Intercept)`)
blup_mun
## mun BLUP
## 1 Chi -3.2150660
## 2 Gig 0.7328851
## 3 HtC 2.8474224
## 4 Jam -0.6709794
## 5 PtR 0.4914815
## 6 RiN -0.5993956
## 7 SnV -2.3350274
## 8 Tam 2.4196258
## 9 ViG 0.8917810
## 10 Yac -0.5627273
#Predichos completos
datos$pred <- predict(modelo_blup)
datos$pred
## [1] 85.60850 87.88766 72.03353 94.13409 85.71846 51.31120 76.39191 86.32342
## [9] 83.87287 86.15202 70.29790 92.39846 83.98283 49.57557 74.65628 84.58779
## [17] 82.99283 85.27199 69.41786 91.51842 83.10279 48.69553 73.77624 83.70775
## [25] 86.94078 89.21994 73.36581 95.46637 87.05074 52.64348 77.72419 87.65570
## [33] 89.05532 91.33447 75.48035 97.58091 89.16528 54.75802 79.83873 89.77024
## [41] 85.53692 87.81607 71.96194 85.64688 51.23962 76.32033 86.69938 88.97853
## [49] 73.12440 86.80934 52.40208 77.48279 87.41430 88.62752 90.90668 75.05255
## [57] 97.15311 88.73748 54.33022 79.41093 89.34244 87.09968 89.37883 73.52470
## [65] 95.62527 87.20964 52.80238 77.88309 87.81459 85.64517 87.92432 72.07020
## [73] 94.17076 85.75513 51.34787 76.42858 86.36009
#Visualizar Blups gen
ggplot(blup_gen, aes(x=reorder(gen, BLUP), y=BLUP)) +
geom_point(size=3) +
geom_hline(yintercept=0, linetype="dashed") +
labs(x = "Genotipo") +
coord_flip()
#Visualizar Blups mun
ggplot(blup_mun, aes(x=reorder(mun, BLUP), y=BLUP)) +
geom_point(size=3) +
geom_hline(yintercept=0, linetype="dashed") +
labs(x = "Municipio") +
coord_flip()
##Componentes de varianza-heredabilidades
vc <- as.data.frame(VarCorr(modelo_blup))
vc
## grp var1 var2 vcov sdcor
## 1 mun (Intercept) <NA> 3.6918174 1.9214103
## 2 gen (Intercept) <NA> 180.8726306 13.4488896
## 3 Residual <NA> <NA> 0.6699132 0.8184823
VarCorr(modelo_blup)
## Groups Name Std.Dev.
## mun (Intercept) 1.92141
## gen (Intercept) 13.44889
## Residual 0.81848
varG <- vc$vcov[vc$grp=="gen"]
varGE <- vc$vcov[vc$grp=="gen:mun"]
varE <- vc$vcov[vc$grp=="Residual"]
e <- 10
r <- 4
# Heredabilidad genotipos y plasticidad
#H2 genotipos
H2g <- varG / (varG + varGE/e + varE/(r*e))
H2g
## numeric(0)
#H2 plasticidad
H2ge <- varGE / (varG + varGE/e + varE/(r*e))
H2ge
## numeric(0)
###ranking genotipos predichos
blup_gen <- ranef(modelo_blup)$gen
blup_gen
## (Intercept)
## CNCH12 5.682404
## CNCH13 7.961560
## FBO1 -7.892568
## FCHI8 14.207994
## FEAR5 5.792363
## FGI4 -28.614893
## FMA7 -3.534182
## FSV1 6.397323
##Predicho de carbono por genotipo
media <- fixef(modelo_blup)[1]
blup_gen$pred <- media + blup_gen[,1]
# Ranking predichos (publicar)
blup_gen[order(-blup_gen$pred),]
## (Intercept) pred
## FCHI8 14.207994 94.73348
## CNCH13 7.961560 88.48705
## FSV1 6.397323 86.92281
## FEAR5 5.792363 86.31785
## CNCH12 5.682404 86.20790
## FMA7 -3.534182 76.99131
## FBO1 -7.892568 72.63292
## FGI4 -28.614893 51.91060
# Visualización ranking predichos
blup_gen$gen <- rownames(blup_gen)
ggplot(blup_gen, aes(x=reorder(gen,pred), y=pred))+
geom_point(size=3)+
coord_flip()+
ylab("Fenotipo predicho (BLUP)")+
xlab("Genotipo")
### Estabilidad con Metan
modelo_metan <- gge(datos, mun, gen, teta)
## Warning: Data imputation used to fill the GxE matrix
modelo_metan
## $teta
## $coordgen
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.5291930 0.2626070 -1.7632742 -0.3667891 1.8993467 -1.4902819
## [2,] 0.7402750 -1.0233096 -0.2803347 -0.0673165 -2.3154278 -1.5947474
## [3,] -0.7354272 -1.5201270 0.9998303 -0.1416089 0.9671007 -0.1157552
## [4,] 1.3214616 1.1259189 0.5873609 2.4262955 0.1335066 0.3428334
## [5,] 0.5402841 1.3176974 1.9603414 -1.8112344 0.1358684 -0.1920216
## [6,] -2.6611458 1.2773039 -0.4009825 0.4828623 -0.6376342 -0.2182060
## [7,] -0.3301461 -1.6826009 0.3993320 0.5314507 0.4512168 1.0718093
## [8,] 0.5955053 0.2425102 -1.5022732 -1.0536596 -0.6339773 2.1963694
## [,7] [,8]
## [1,] -0.4027506 1.170064
## [2,] -0.0251020 1.170064
## [3,] 2.1825004 1.170064
## [4,] 0.4495869 1.170064
## [5,] -0.6134081 1.170064
## [6,] -0.1483500 1.170064
## [7,] -2.2020272 1.170064
## [8,] 0.7595507 1.170064
##
## $coordenv
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 38.11635 2.00349663 -0.5102912 -0.7966537 0.5131806 0.091670554
## [2,] 34.29211 0.04167356 0.2004954 -0.2562565 0.4252528 -0.365269375
## [3,] 34.13375 -0.52671105 0.3639158 0.9876035 0.4043181 0.213004594
## [4,] 35.88459 0.05176954 0.6272620 -0.4386836 -0.3104298 -0.006472464
## [5,] 35.40430 -0.43343584 0.6583729 0.2400356 -0.2003434 -0.059676859
## [6,] 36.80466 2.42834546 -0.3824267 0.6231082 -0.5602712 0.012790250
## [7,] 35.20830 -1.60349012 -1.9634013 0.1494122 -0.1550262 -0.039992456
## [8,] 34.46362 -0.44105653 0.2330611 0.6955156 0.4333166 0.015194207
## [9,] 35.33028 -0.65168705 0.5837500 0.1147031 -0.4000101 -0.063680469
## [10,] 36.21328 -1.13279566 0.2452658 -1.2020486 -0.1224414 0.186938917
## [,7] [,8]
## [1,] 0.039945164 -1.140879e-15
## [2,] -0.066740140 2.917442e-15
## [3,] -0.110064762 9.041787e-16
## [4,] -0.049065450 -1.159635e-14
## [5,] -0.002916231 1.841941e-15
## [6,] -0.018113561 3.311919e-15
## [7,] -0.009913124 -1.851959e-15
## [8,] 0.135267121 -2.431929e-15
## [9,] 0.087197877 2.305071e-15
## [10,] -0.009385303 5.776335e-15
##
## $eigenvalues
## [1] 1.125907e+02 3.854115e+00 2.390894e+00 2.068370e+00 1.202028e+00
## [6] 4.814505e-01 2.167837e-01 1.440454e-14
##
## $totalvar
## [1] 12703.24
##
## $varexpl
## [1] 99.79 0.12 0.04 0.03 0.01 0.00 0.00 0.00
##
## $labelgen
## [1] "CNCH12" "CNCH13" "FBO1" "FCHI8" "FEAR5" "FGI4" "FMA7" "FSV1"
##
## $labelenv
## [1] "Chi" "Gig" "HtC" "Jam" "PtR" "RiN" "SnV" "Tam" "ViG" "Yac"
##
## $labelaxes
## [1] "PC1" "PC2" "PC3" "PC4" "PC5" "PC6" "PC7" "PC8"
##
## $ge_mat
## Chi Gig HtC Jam PtR RiN SnV Tam
## CNCH12 6.8625 5.825 5.2625 5.287339 5.161762 5.8875 6.4625 5.5
## CNCH13 7.5625 7.525 7.3625 8.187339 8.161762 7.8875 8.6625 7.5
## FBO1 -9.3375 -7.475 -7.2375 -7.912661 -7.538238 -9.6125 -7.7375 -7.2
## FCHI8 15.2625 13.525 14.2625 14.116203 14.267669 15.8875 13.2625 14.2
## FEAR5 7.1625 5.925 5.0625 6.487339 5.861762 6.3875 3.8625 5.2
## FGI4 -30.0375 -27.675 -27.6375 -28.912661 -28.638238 -28.4125 -28.6375 -27.9
## FMA7 -4.9375 -3.475 -2.7375 -3.612661 -3.238238 -4.9125 -2.9375 -3.1
## FSV1 7.4625 5.825 5.6625 6.359760 5.961762 6.8875 7.0625 5.8
## ViG Yac
## CNCH12 5.0625 5.55
## CNCH13 8.3625 8.45
## FBO1 -7.4375 -7.45
## FCHI8 14.0625 13.25
## FEAR5 5.7625 6.25
## FGI4 -28.6375 -29.75
## FMA7 -3.2375 -3.15
## FSV1 6.0625 6.85
##
## $centering
## [1] "environment"
##
## $scaling
## [1] "none"
##
## $svp
## [1] "environment"
##
## $d
## [1] 0.3021659
##
## $grand_mean
## [1] 80.52384
##
## $mean_gen
## CNCH12 CNCH13 FBO1 FCHI8 FEAR5 FGI4 FMA7 FSV1
## 86.21000 88.49000 72.63000 94.73348 86.32000 51.90000 76.99000 86.91724
##
## $mean_env
## Chi Gig HtC Jam PtR RiN SnV Tam
## 77.23750 81.27500 83.43750 79.81266 81.03824 79.91250 78.13750 83.00000
## ViG Yac
## 81.43750 79.95000
##
## $scale_val
## Chi Gig HtC Jam PtR RiN SnV Tam
## 14.43230 12.96355 12.91023 13.56672 13.38543 13.94545 13.34219 13.03117
## ViG Yac
## 13.35867 13.70214
##
## attr(,"class")
## [1] "gge"
##
## attr(,"class")
## [1] "gge"
#Grafica metan
plot(modelo_metan)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the metan package.
## Please report the issue at <https://github.com/nepem-ufsc/metan/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the metan package.
## Please report the issue at <https://github.com/nepem-ufsc/metan/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
### Selección ideotípica
##Integración de variables
##Blup gen
blup_gen <- ranef(modelo_blup)$gen
media <- fixef(modelo_blup)[1]
blup_gen$BLUP_C <- media + blup_gen[,1]
blup_gen$gen <- rownames(blup_gen)
blup_gen <- blup_gen[,c("gen","BLUP_C")]
blup_gen
## gen BLUP_C
## CNCH12 CNCH12 86.20790
## CNCH13 CNCH13 88.48705
## FBO1 FBO1 72.63292
## FCHI8 FCHI8 94.73348
## FEAR5 FEAR5 86.31785
## FGI4 FGI4 51.91060
## FMA7 FMA7 76.99131
## FSV1 FSV1 86.92281
##Plasticidad usando Fisher environments (joint regression)
#índice (creando valores de x para definir env = promedio de tasas en c/parcela)
indice_env <- datos %>%
group_by(mun) %>%
summarise(env = mean(teta))
datos <- left_join(datos, indice_env, by="mun")
#visualización Normas de reacción joint regression env
ggplot(datos, aes(x = env, y = teta,
color = gen)) +
geom_smooth(method = "lm", se = FALSE) +
labs(x = "Ambiente (local)",
y = expression(Fenotipo)) +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
#visualización Normas de reacción clima local
ggplot(datos, aes(x = E, y = teta,
color = gen)) +
geom_smooth(method = "lm", se = FALSE) +
labs(x = "Ambiente (Estrés)",
y = expression(Fenotipo)) +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
## plasticidad joint
# modelo factores fijos
mod_plas_lm <- lm(teta ~ gen*env,
data=datos)
y=emtrends(mod_plas_lm, "gen", var = "env")
multcomp::cld(y, Letters = LETTERS)
## gen env.trend SE df lower.CL upper.CL .group
## FEAR5 0.572 0.174 61 0.224 0.920 A
## CNCH12 0.603 0.174 61 0.255 0.951 A
## CNCH13 0.619 0.174 61 0.272 0.967 A
## FSV1 0.656 0.206 61 0.244 1.068 A
## FMA7 0.863 0.174 61 0.515 1.210 A
## FBO1 0.896 0.174 61 0.548 1.244 A
## FCHI8 0.934 0.213 61 0.509 1.359 A
## FGI4 0.951 0.174 61 0.603 1.299 A
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 8 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.
# modelo blup factores aleatorios
modelo_plasticidad <- lmer(teta ~ env +
(env|gen) +
(1|mun),
data=datos)
## boundary (singular) fit: see help('isSingular')
pend <- ranef(modelo_plasticidad)$gen
pend$gen <- rownames(pend)
plasticidad <- pend[,c("gen","env")]
colnames(plasticidad)[2] <- "Pendiente"
#plasticidad Estrés
# modelo factores fijos
mod_plas2_lm <- lm(teta ~ gen*E,
data=datos)
z<-emtrends(mod_plas2_lm, "gen", var = "E")
multcomp::cld(z, Letters = LETTERS)
## gen E.trend SE df lower.CL upper.CL .group
## FSV1 17.1 8.49 61 0.16 34.1 A
## CNCH13 19.2 8.38 61 2.50 36.0 A
## FEAR5 20.0 8.38 61 3.25 36.8 A
## CNCH12 21.7 8.38 61 4.98 38.5 A
## FMA7 27.5 8.38 61 10.72 44.2 A
## FBO1 28.6 8.38 61 11.82 45.3 A
## FCHI8 30.0 8.86 61 12.27 47.7 A
## FGI4 35.6 8.38 61 18.82 52.3 A
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 8 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.