Entrega Final

Curso rápido de RLEM

Dados

A seguir, criamos o pedrigree conforme solicitado, com 20 indivíduos. Desses, 12 foram fenotipados e os demais não o foram. A tabela a seguir expõe os dados criados:

Code

ped <- data.frame(
  ped = 1,
  id = LETTERS[1:20],
  sex = c(2,1,2,1,
          1,2,1,2,2,2,1,2,2,
          1,2,1,1,1,2,1),
  mother = c(rep(NA,6),rep("A", 4), rep("C", 3), rep("F",2), rep("J", 4), "F"),
  father = c(rep(NA, 6), rep("B", 4), rep("D", 3), rep("G", 2), rep("K", 4), "E"),
  pheno = c(0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,1,0,1,0,1),
  rep = "A",
  y = c(rep(NA, 4), 3, 6.8, 8, 9.1, 7.2, 5, 4.1, 5, 3.6, NA, NA, 7, NA, 6.5, NA, 6)
)

ped  %>%
  knitr::kable(
    align = "c",
    booktabs = TRUE,
    longtable = TRUE,
    linesep = "",
    escape = FALSE
    )  %>%
  kableExtra::kable_styling(
      position = "center",
      latex_options = c("striped", "repeat_header"),
      stripe_color = "gray!15")

ped	id	sex	mother	father	pheno	rep	y
1	A	2	NA	NA	0	A	NA
1	B	1	NA	NA	0	A	NA
1	C	2	NA	NA	0	A	NA
1	D	1	NA	NA	0	A	NA
1	E	1	NA	NA	1	A	3.0
1	F	2	NA	NA	1	A	6.8
1	G	1	A	B	1	A	8.0
1	H	2	A	B	1	A	9.1
1	I	2	A	B	1	A	7.2
1	J	2	A	B	1	A	5.0
1	K	1	C	D	1	A	4.1
1	L	2	C	D	1	A	5.0
1	M	2	C	D	1	A	3.6
1	N	1	F	G	0	A	NA
1	O	2	F	G	0	A	NA
1	P	1	J	K	1	A	7.0
1	Q	1	J	K	0	A	NA
1	R	1	J	K	1	A	6.5
1	S	2	J	K	0	A	NA
1	T	1	F	E	1	A	6.0

Representamos a seguir o pedigree criado:

Code

ped_plot <- pedigree(id = ped$id,
                     dadid = ped$father,
                     momid = ped$mother,
                     sex = ped$sex,
                     famid = ped$ped,
                     affected = ped$pheno)


plot(ped_plot['1'], col=ifelse(ped$pheno==0,"black", "green"), cex = .7)

RLEM

O próximo passo para a geração do modelo de RLEM é a criação da matriz aditiva, uma matriz de covariâncias advindas de parentesco. A matriz é criada no bloco de código a seguir:

Code

A <- ped %>%
  dplyr::select("ID" = id, "Dam" = mother, "Sire" = father) %>%
  mutate(across(.cols = Dam:Sire, ~if_else(is.na(.), "0", .))
  )%>%
  Amatrix(ploidy = 2)

Nosso modelo, considerando apenas os individuos fenotipados, foi construído utilizando sexo como efeito fixo e os indivíduos como efeitos aleatórios. O resumo do modelo e BLUPs são apresentados a seguir.

Code

phe <- ped %>%
  dplyr::select("ID" = id, sex, y) %>%
  mutate(sex = as.factor(sex))
  
fitI <- mmer(y~sex, # efeito fixo
     random= ~ ID,
     rcov= ~ units, # todas as unidades fazem parte do erro
     data=phe)

Code

summary(fitI) #resumao do fitI

============================================================
         Multivariate Linear Mixed Model fit by REML         
**********************  sommer 4.2  ********************** 
============================================================
         logLik      AIC      BIC Method Converge
Value -5.427003 14.85401 15.82382     NR     TRUE
============================================================
Variance-Covariance components:
          VarComp VarCompSE Zratio Constraint
ID.y-y      0.615    0.8251 0.7454   Positive
units.y-y   3.075    0.8251 3.7268   Positive
============================================================
Fixed effects:
  Trait      Effect Estimate Std.Error t.value
1     y (Intercept)    5.767    0.7842  7.3532
2     y        sex2    0.350    1.1091  0.3156
============================================================
Groups and observations:
    y
ID 12
============================================================
Use the '$' sign to access results and parameters

Code

fitI$U$ID$y #BLUPs

        IDE         IDF         IDG         IDH         IDI         IDJ 
-0.46111111  0.11388889  0.37222222  0.49722222  0.18055556 -0.18611111 
        IDK         IDL         IDM         IDP         IDR         IDT 
-0.27777778 -0.18611111 -0.41944444  0.20555556  0.12222222  0.03888889

De fato, a prole do casal AB parece contribuir para maiores níveis da variável \(y\) medida.

Consideremos agora o modelo para toda a família, ponderando o modelo anterior pela matriz de covariâncias de pedigree. Os efeitos fixos e os aleatórios são os mesmos.

Code

#Modelo conhecendo o pedigree (matriz Aditiva)
fitA <- mmer(y~sex,
             random= ~ vsr(ID, Gu=A),
             rcov= ~ units,
             data=phe)

Code

summary(fitA) #resumao do fitA

============================================================
         Multivariate Linear Mixed Model fit by REML         
**********************  sommer 4.2  ********************** 
============================================================
         logLik      AIC      BIC Method Converge
Value -4.692781 13.38556 14.35538     NR     TRUE
============================================================
Variance-Covariance components:
          VarComp VarCompSE Zratio Constraint
u:ID.y-y    2.583     3.425 0.7541   Positive
units.y-y   1.324     2.152 0.6152   Positive
============================================================
Fixed effects:
  Trait      Effect Estimate Std.Error t.value
1     y (Intercept)   5.8700    0.9809   5.984
2     y        sex2  -0.2907    1.0095  -0.288
============================================================
Groups and observations:
      y
u:ID 20
============================================================
Use the '$' sign to access results and parameters

Por fim, temos efeitos fixos estimados \(\beta_{0M} =\) 5.870002 e \(\beta_{0F} =\) -0.2907265 e os efeitos aleatórios, bem como os valores \(\hat{y}\), são apresentados na tabela a seguir:

Code

#estimacao dos valores
tiny <- ped %>% dplyr::select(id, sex, y)

yhat <- data.frame(individuo = names(fitA$U$`u:ID`$y),
         BLUP = fitA$U$`u:ID`$y) %>%
  arrange(individuo) %>%
  inner_join(tiny, by = c("individuo" = "id"))%>%
  mutate(estimacao = if_else(sex == 1, BLUP + fitA$Beta$Estimate[1], BLUP + fitA$Beta$Estimate[2])
         )



row.names(yhat) <- NULL

yhat%>%
  knitr::kable(
    align = "c",
    booktabs = TRUE,
    longtable = TRUE,
    linesep = "",
    escape = TRUE,
    digits = 4
    ) %>%
  kableExtra::kable_styling(
      position = "center",
      latex_options = c("striped", "repeat_header"),
      stripe_color = "gray!15")

individuo	BLUP	sex	y	estimacao
A	1.2012	2	NA	0.9104
B	1.2012	1	NA	7.0712
C	-0.7528	2	NA	-1.0435
D	-0.7528	1	NA	5.1172
E	-1.8005	1	3.0	4.0695
F	0.9038	2	6.8	0.6130
G	1.6597	1	8.0	7.5297
H	2.3464	2	9.1	2.0556
I	1.4083	2	7.2	1.1176
J	0.5914	2	5.0	0.3007
K	-0.9858	1	4.1	4.8842
L	-0.6671	2	5.0	-0.9579
M	-1.3583	2	3.6	-1.6491
N	1.2818	1	NA	7.1518
O	1.2818	2	NA	0.9910
P	0.4581	1	7.0	6.3281
Q	-0.1972	1	NA	5.6728
R	0.2112	1	6.5	6.0812
S	-0.1972	2	NA	-0.4879
T	-0.1628	1	6.0	5.7072

Finalmente, comparamos os valores reais com os estimados, auxiliados pela reta de regressão ajustada.

Code

yhat %>%
  dplyr::filter(!is.na(y) == T) %>%
  ggplot(aes(y, estimacao))+
  geom_point()+
  geom_smooth(method = "lm",se = FALSE)+
  theme_bw() +
  theme(axis.title.y=element_text(colour="black", size=12),
        axis.title.x = element_text(colour="black", size=12),
        axis.text = element_text(colour = '#383A46', size=9.5),
        panel.border = element_blank(),
        axis.line = element_line(colour = '#383A46'),
        panel.grid.minor = element_blank())

`geom_smooth()` using formula = 'y ~ x'

ped	id	sex	mother	father	pheno	rep	y
1	A	2	NA	NA	0	A	NA
1	B	1	NA	NA	0	A	NA
1	C	2	NA	NA	0	A	NA
1	D	1	NA	NA	0	A	NA
1	E	1	NA	NA	1	A	3.0
1	F	2	NA	NA	1	A	6.8
1	G	1	A	B	1	A	8.0
1	H	2	A	B	1	A	9.1
1	I	2	A	B	1	A	7.2
1	J	2	A	B	1	A	5.0
1	K	1	C	D	1	A	4.1
1	L	2	C	D	1	A	5.0
1	M	2	C	D	1	A	3.6
1	N	1	F	G	0	A	NA
1	O	2	F	G	0	A	NA
1	P	1	J	K	1	A	7.0
1	Q	1	J	K	0	A	NA
1	R	1	J	K	1	A	6.5
1	S	2	J	K	0	A	NA
1	T	1	F	E	1	A	6.0

ped	id	sex	mother	father	pheno	rep	y
1	A	2	NA	NA	0	A	NA
1	B	1	NA	NA	0	A	NA
1	C	2	NA	NA	0	A	NA
1	D	1	NA	NA	0	A	NA
1	E	1	NA	NA	1	A	3.0
1	F	2	NA	NA	1	A	6.8
1	G	1	A	B	1	A	8.0
1	H	2	A	B	1	A	9.1
1	I	2	A	B	1	A	7.2
1	J	2	A	B	1	A	5.0
1	K	1	C	D	1	A	4.1
1	L	2	C	D	1	A	5.0
1	M	2	C	D	1	A	3.6
1	N	1	F	G	0	A	NA
1	O	2	F	G	0	A	NA
1	P	1	J	K	1	A	7.0
1	Q	1	J	K	0	A	NA
1	R	1	J	K	1	A	6.5
1	S	2	J	K	0	A	NA
1	T	1	F	E	1	A	6.0

ped	id	sex	mother	father	pheno	rep	y
1	A	2	NA	NA	0	A	NA
1	B	1	NA	NA	0	A	NA
1	C	2	NA	NA	0	A	NA
1	D	1	NA	NA	0	A	NA
1	E	1	NA	NA	1	A	3.0
1	F	2	NA	NA	1	A	6.8
1	G	1	A	B	1	A	8.0
1	H	2	A	B	1	A	9.1
1	I	2	A	B	1	A	7.2
1	J	2	A	B	1	A	5.0
1	K	1	C	D	1	A	4.1
1	L	2	C	D	1	A	5.0
1	M	2	C	D	1	A	3.6
1	N	1	F	G	0	A	NA
1	O	2	F	G	0	A	NA
1	P	1	J	K	1	A	7.0
1	Q	1	J	K	0	A	NA
1	R	1	J	K	1	A	6.5
1	S	2	J	K	0	A	NA
1	T	1	F	E	1	A	6.0