1 Exercício 1

Dados:

Os pais forma codificados de 1 a 5, e o filhos de 11 a 16, sendo s = macho e d = fêmea.

dados
##   Animal s d Sexo PNascer  P42 Ovos
## 1     11 0 0    0       0    0    0
## 2     12 1 0    2      52 2410  205
## 3     13 1 2    1      54 2450    0
## 4     14 1 2    2      53 2410  215
## 5     15 3 4    1      55 2520    0
## 6     16 5 4    1      51 2215    0

1.1 Modelos:

Modelo 1: Característica 1= Peso aos 42 dias

\[ Y{ij} = \mu + S_i + b(P_{ij} - \overline{P}) + a_{ij} + e_{ij} \]

Modelo 2: Característica 2 = Número de ovos

\[ Y{ij} = \mu + a_i + e_{ij} \]

1.2 Obter a matriz dos numeradores dos coeficientes de parentesco (A) usando a covariância entre os indivíduos;

install.packages("AGHmatrix", repos = "http://cran.us.r-project.org")
## package 'AGHmatrix' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\ferre\AppData\Local\Temp\Rtmp4KDwsp\downloaded_packages
library("AGHmatrix")
Acov <- Amatrix(P, ploidy = 2)
## Verifying conflicting data... 
## Organizing data... 
## Your data was chronologically organized with success. 
## Constructing matrix A using ploidy = 2 
## Completed! Time = 0.0003333333  minutes
Acov
##      1   2   3   4   5 11   12   13   14   15   16
## 1  1.0 0.0 0.0 0.0 0.0  0 0.50 0.50 0.50 0.00 0.00
## 2  0.0 1.0 0.0 0.0 0.0  0 0.00 0.50 0.50 0.00 0.00
## 3  0.0 0.0 1.0 0.0 0.0  0 0.00 0.00 0.00 0.50 0.00
## 4  0.0 0.0 0.0 1.0 0.0  0 0.00 0.00 0.00 0.50 0.50
## 5  0.0 0.0 0.0 0.0 1.0  0 0.00 0.00 0.00 0.00 0.50
## 11 0.0 0.0 0.0 0.0 0.0  1 0.00 0.00 0.00 0.00 0.00
## 12 0.5 0.0 0.0 0.0 0.0  0 1.00 0.25 0.25 0.00 0.00
## 13 0.5 0.5 0.0 0.0 0.0  0 0.25 1.00 0.50 0.00 0.00
## 14 0.5 0.5 0.0 0.0 0.0  0 0.25 0.50 1.00 0.00 0.00
## 15 0.0 0.0 0.5 0.5 0.0  0 0.00 0.00 0.00 1.00 0.25
## 16 0.0 0.0 0.0 0.5 0.5  0 0.00 0.00 0.00 0.25 1.00

Matriz A

1.3 Obter a matriz dos numeradores dos coeficientes de parentesco (A) e sua inversa (A-¹) usando as regras de Quaas (Henderson);

Literatura: Henderson, C. R. 1976. Simple Method for Computing the Inverse of a Numerator Relationship Matrix Used in Prediction of Breeding Values. Biometrics 32:69-83. Function: http://morotalab.org/Mrode2005/relmat/relmat.html

Matriz A e A-¹ (Ai)

s = ( c (0 , 1, 1, 1, 3, 5))
d = (c ( 0, 0, 2, 2, 4, 4))
Aquass <- createA(s,d)
Aquass
##      [,1] [,2]  [,3]  [,4]   [,5]   [,6]
## [1,] 1.00 0.50 0.750 0.750 0.7500 0.7500
## [2,] 0.50 1.00 0.750 0.750 0.7500 0.7500
## [3,] 0.75 0.75 1.250 0.750 1.0000 0.8750
## [4,] 0.75 0.75 0.750 1.250 1.0000 1.1250
## [5,] 0.75 0.75 1.000 1.000 1.3750 1.1875
## [6,] 0.75 0.75 0.875 1.125 1.1875 1.5000
Aquass_inv <- solve(Aquass)
Aquass_inv
##            [,1]       [,2]       [,3]       [,4]       [,5]          [,6]
## [1,]  2.3333333  0.3333333 -1.0000000 -1.0000000  0.0000000  5.551115e-17
## [2,]  0.3333333  2.3333333 -1.0000000 -1.0000000  0.0000000  5.551115e-17
## [3,] -1.0000000 -1.0000000  2.6666667  0.6666667 -1.3333333 -1.110223e-16
## [4,] -1.0000000 -1.0000000  0.6666667  3.3939394 -0.6060606 -1.454545e+00
## [5,]  0.0000000  0.0000000 -1.3333333 -0.6060606  3.3939394 -1.454545e+00
## [6,]  0.0000000  0.0000000  0.0000000 -1.4545455 -1.4545455  2.909091e+00

1.4 Montar o vetor y; as matrizes X, Z, R e G, o sistema de equações de modelos mistos e obter as soluções para as equações de modelos mistos usando análise de características múltiplas para peso aos 42 dias e número de ovos.

Usar característica dentro de animal:

\[ G= A\otimes G_0 \]

Reorganizar o banco de dados para a análise multi características.

Passos:

  1. Instalar o pacotesreshape2
  2. Fundir duas colunas em linhas O comando melt realiza basicamente o agrupamento de todas as colunas não listadas, conforme o seguinte (neste caso agrupamos por “P42” e “Ovos)”:
  3. Ordenar o dados_melt por animal
dados_melt <- reshape2::melt (dados, id.vars=c("Animal", "s", "d", "Sexo", "PNascer"))
dadosc<- dados_melt[order(dados_melt$Animal),]
dadosc
##    Animal s d Sexo PNascer variable value
## 1      11 0 0    0       0      P42     0
## 7      11 0 0    0       0     Ovos     0
## 2      12 1 0    2      52      P42  2410
## 8      12 1 0    2      52     Ovos   205
## 3      13 1 2    1      54      P42  2450
## 9      13 1 2    1      54     Ovos     0
## 4      14 1 2    2      53      P42  2410
## 10     14 1 2    2      53     Ovos   215
## 5      15 3 4    1      55      P42  2520
## 11     15 3 4    1      55     Ovos     0
## 6      16 5 4    1      51      P42  2215
## 12     16 5 4    1      51     Ovos     0

Passo 2: Escrever o vetor Y

Yc <-as.matrix(dadosc$value)
Yc
##       [,1]
##  [1,]    0
##  [2,]    0
##  [3,] 2410
##  [4,]  205
##  [5,] 2450
##  [6,]    0
##  [7,] 2410
##  [8,]  215
##  [9,] 2520
## [10,]    0
## [11,] 2215
## [12,]    0

Passo 3: Escrever a matriz X

  1. Criar a coluna para a média da característica 1 (m1), para a média da característica 2(m2), para a codificação do sexo (-1, 0 e 1) de acordo com os modelos (s1)e para a covariável (p1).
m1<- dadosc$value
m1[dadosc$variable == "P42" & dadosc$value > 0] = 1
m1[dadosc$variable == "Ovos"] = 0
m2<- dadosc$value
m2[dadosc$variable == "Ovos" & dadosc$value > 0] = 1
m2[dadosc$variable == "P42"] = 0
s1 <- dadosc$Sexo
s1 [dadosc$Sexo == 2 & dadosc$variable == "P42"] = -1
s1 [dadosc$Sexo == 1 & dadosc$variable == "P42"] = 1
s1 [dadosc$Sexo == 2 & dadosc$variable == "Ovos"] = 0
s1 [dadosc$Sexo == 1 & dadosc$variable == "Ovos"] = 0
soma <-sum (dadosc$PNascer)
n <- length(dadosc$PNascer[dadosc$PNascer !=0])
md <- soma/n
p1 <- dadosc$PNascer - md
p1 [p1 < -15]<- 0
Xc <- data.frame(m1,m2,s1,p1)
Xc <-as.matrix(Xc)
Xc
##       m1 m2 s1 p1
##  [1,]  0  0  0  0
##  [2,]  0  0  0  0
##  [3,]  1  0 -1 -1
##  [4,]  0  1  0 -1
##  [5,]  1  0  1  1
##  [6,]  0  0  0  1
##  [7,]  1  0 -1  0
##  [8,]  0  1  0  0
##  [9,]  1  0  1  2
## [10,]  0  0  0  2
## [11,]  1  0  1 -2
## [12,]  0  0  0 -2

Passo 4: Escrever o vetor Z

Zc <- read.csv("C:\\Users\\ferre\\Dropbox\\ZOO 760\\Relatorio_2\\matriz_Z_lista2.txt",h=F,sep=" ")
Zc <- as.matrix(Zc)
Zc
##       V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
##  [1,]  0  0  0  0  0  0  0  0  0   0   0   0
##  [2,]  0  0  0  0  0  0  0  0  0   0   0   0
##  [3,]  0  0  1  0  0  0  0  0  0   0   0   0
##  [4,]  0  0  0  1  0  0  0  0  0   0   0   0
##  [5,]  0  0  0  0  1  0  0  0  0   0   0   0
##  [6,]  0  0  0  0  0  0  0  0  0   0   0   0
##  [7,]  0  0  0  0  0  0  1  0  0   0   0   0
##  [8,]  0  0  0  0  0  0  0  1  0   0   0   0
##  [9,]  0  0  0  0  0  0  0  0  1   0   0   0
## [10,]  0  0  0  0  0  0  0  0  0   0   0   0
## [11,]  0  0  0  0  0  0  0  0  0   0   1   0
## [12,]  0  0  0  0  0  0  0  0  0   0   0   0

Passo 5: Obter a matriz R e G e suas respectivas inversas

Gzero <- rbind(c (10000, -100), c (-100, 80))
Gzero <-as.matrix(Gzero)
Gc <- Aquass %x% Gzero
Gc
##        [,1] [,2]  [,3] [,4]    [,5]   [,6]    [,7]   [,8]     [,9]   [,10]
##  [1,] 10000 -100  5000  -50  7500.0  -75.0  7500.0  -75.0  7500.00  -75.00
##  [2,]  -100   80   -50   40   -75.0   60.0   -75.0   60.0   -75.00   60.00
##  [3,]  5000  -50 10000 -100  7500.0  -75.0  7500.0  -75.0  7500.00  -75.00
##  [4,]   -50   40  -100   80   -75.0   60.0   -75.0   60.0   -75.00   60.00
##  [5,]  7500  -75  7500  -75 12500.0 -125.0  7500.0  -75.0 10000.00 -100.00
##  [6,]   -75   60   -75   60  -125.0  100.0   -75.0   60.0  -100.00   80.00
##  [7,]  7500  -75  7500  -75  7500.0  -75.0 12500.0 -125.0 10000.00 -100.00
##  [8,]   -75   60   -75   60   -75.0   60.0  -125.0  100.0  -100.00   80.00
##  [9,]  7500  -75  7500  -75 10000.0 -100.0 10000.0 -100.0 13750.00 -137.50
## [10,]   -75   60   -75   60  -100.0   80.0  -100.0   80.0  -137.50  110.00
## [11,]  7500  -75  7500  -75  8750.0  -87.5 11250.0 -112.5 11875.00 -118.75
## [12,]   -75   60   -75   60   -87.5   70.0  -112.5   90.0  -118.75   95.00
##          [,11]   [,12]
##  [1,]  7500.00  -75.00
##  [2,]   -75.00   60.00
##  [3,]  7500.00  -75.00
##  [4,]   -75.00   60.00
##  [5,]  8750.00  -87.50
##  [6,]   -87.50   70.00
##  [7,] 11250.00 -112.50
##  [8,]  -112.50   90.00
##  [9,] 11875.00 -118.75
## [10,]  -118.75   95.00
## [11,] 15000.00 -150.00
## [12,]  -150.00  120.00
Gic <- solve(Gc)
Gic <- as.matrix(Gic)
Ic <- as.matrix (diag(6))
Rzero <- rbind(c (10000, -250), c (-250, 200))
Rc <- Ic %x% Rzero
Rc
##        [,1] [,2]  [,3] [,4]  [,5] [,6]  [,7] [,8]  [,9] [,10] [,11] [,12]
##  [1,] 10000 -250     0    0     0    0     0    0     0     0     0     0
##  [2,]  -250  200     0    0     0    0     0    0     0     0     0     0
##  [3,]     0    0 10000 -250     0    0     0    0     0     0     0     0
##  [4,]     0    0  -250  200     0    0     0    0     0     0     0     0
##  [5,]     0    0     0    0 10000 -250     0    0     0     0     0     0
##  [6,]     0    0     0    0  -250  200     0    0     0     0     0     0
##  [7,]     0    0     0    0     0    0 10000 -250     0     0     0     0
##  [8,]     0    0     0    0     0    0  -250  200     0     0     0     0
##  [9,]     0    0     0    0     0    0     0    0 10000  -250     0     0
## [10,]     0    0     0    0     0    0     0    0  -250   200     0     0
## [11,]     0    0     0    0     0    0     0    0     0     0 10000  -250
## [12,]     0    0     0    0     0    0     0    0     0     0  -250   200
Ric <- solve(Rc)
Ric <- as.matrix(Ric)

Passo 6: Solucionar as Equações de Modelos Mistos (EMM)

um <- t(Xc) %*% Ric %*% Xc
dois<- t(Zc) %*% Ric %*% Xc
tres<- t(Xc) %*% Ric %*% Zc
quatro<- (t(Zc) %*% Ric %*% Zc) + Gic
cinco<- t(Xc) %*% Ric %*% Yc
seis<- t(Zc) %*% Ric %*% Yc

LHSc <-cbind(rbind(um,dois),rbind(tres, quatro))
RHSc <-rbind(cinco, seis)
paramc <-MASS::ginv(LHSc) %*% RHSc
paramc
##               [,1]
##  [1,] 2408.2677585
##  [2,]  211.1301764
##  [3,]  -13.2972820
##  [4,]    2.7899189
##  [5,]    2.2504639
##  [6,]    0.3036598
##  [7,]   -2.2504639
##  [8,]   -0.3036598
##  [9,]   22.5906260
## [10,]   -0.2259063
## [11,]  -18.0896981
## [12,]    0.8332258
## [13,]   15.8822765
## [14,]    0.1673416
## [15,]  -44.6616497
## [16,]    0.9358631

Refinando…

paramc<- as.data.frame(paramc)
paramc<- tibble::rownames_to_column(paramc, var = "Parametros")
paramc$Parametros <- c("Média P42", "Média Ovos", " Sexo", "Covariável", "A1", "A1", "A2", "A2",
                    "A3", "A3", "A4", "A4", "A5", "A5", "A6", "A6")
colnames(paramc)[2] <- "Efeitos"
paramc
##    Parametros      Efeitos
## 1   Média P42 2408.2677585
## 2  Média Ovos  211.1301764
## 3        Sexo  -13.2972820
## 4  Covariável    2.7899189
## 5          A1    2.2504639
## 6          A1    0.3036598
## 7          A2   -2.2504639
## 8          A2   -0.3036598
## 9          A3   22.5906260
## 10         A3   -0.2259063
## 11         A4  -18.0896981
## 12         A4    0.8332258
## 13         A5   15.8822765
## 14         A5    0.1673416
## 15         A6  -44.6616497
## 16         A6    0.9358631

1.5 Montar o vetor y; as matrizes X, Z, R e G, o sistema de equações de modelos mistos e obter as soluções para as equações de modelos mistos usando análise de características múltiplas para peso aos 42 dias e número de ovos.

Usar animal dentro de característica:

\[ G= G_0 \otimes A \] Passo 1: Obter a matriz G;

Gd <- Gzero %x% Aquass 
Gd
##        [,1]  [,2]    [,3]    [,4]     [,5]     [,6] [,7] [,8]   [,9]  [,10]
##  [1,] 10000  5000  7500.0  7500.0  7500.00  7500.00 -100  -50  -75.0  -75.0
##  [2,]  5000 10000  7500.0  7500.0  7500.00  7500.00  -50 -100  -75.0  -75.0
##  [3,]  7500  7500 12500.0  7500.0 10000.00  8750.00  -75  -75 -125.0  -75.0
##  [4,]  7500  7500  7500.0 12500.0 10000.00 11250.00  -75  -75  -75.0 -125.0
##  [5,]  7500  7500 10000.0 10000.0 13750.00 11875.00  -75  -75 -100.0 -100.0
##  [6,]  7500  7500  8750.0 11250.0 11875.00 15000.00  -75  -75  -87.5 -112.5
##  [7,]  -100   -50   -75.0   -75.0   -75.00   -75.00   80   40   60.0   60.0
##  [8,]   -50  -100   -75.0   -75.0   -75.00   -75.00   40   80   60.0   60.0
##  [9,]   -75   -75  -125.0   -75.0  -100.00   -87.50   60   60  100.0   60.0
## [10,]   -75   -75   -75.0  -125.0  -100.00  -112.50   60   60   60.0  100.0
## [11,]   -75   -75  -100.0  -100.0  -137.50  -118.75   60   60   80.0   80.0
## [12,]   -75   -75   -87.5  -112.5  -118.75  -150.00   60   60   70.0   90.0
##         [,11]   [,12]
##  [1,]  -75.00  -75.00
##  [2,]  -75.00  -75.00
##  [3,] -100.00  -87.50
##  [4,] -100.00 -112.50
##  [5,] -137.50 -118.75
##  [6,] -118.75 -150.00
##  [7,]   60.00   60.00
##  [8,]   60.00   60.00
##  [9,]   80.00   70.00
## [10,]   80.00   90.00
## [11,]  110.00   95.00
## [12,]   95.00  120.00
Gid <- solve(Gd)
Gid <- as.matrix(Gid)

Passo 2: Obter a matriz R;

Id <- as.matrix (diag(6))
Rd <- Rzero %x% Id
Rd
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6] [,7] [,8] [,9] [,10] [,11] [,12]
##  [1,] 10000     0     0     0     0     0 -250    0    0     0     0     0
##  [2,]     0 10000     0     0     0     0    0 -250    0     0     0     0
##  [3,]     0     0 10000     0     0     0    0    0 -250     0     0     0
##  [4,]     0     0     0 10000     0     0    0    0    0  -250     0     0
##  [5,]     0     0     0     0 10000     0    0    0    0     0  -250     0
##  [6,]     0     0     0     0     0 10000    0    0    0     0     0  -250
##  [7,]  -250     0     0     0     0     0  200    0    0     0     0     0
##  [8,]     0  -250     0     0     0     0    0  200    0     0     0     0
##  [9,]     0     0  -250     0     0     0    0    0  200     0     0     0
## [10,]     0     0     0  -250     0     0    0    0    0   200     0     0
## [11,]     0     0     0     0  -250     0    0    0    0     0   200     0
## [12,]     0     0     0     0     0  -250    0    0    0     0     0   200
Rid <- solve(Rd)
Rid <- as.matrix(Rid)

Passo 3: Solucionar as Equações de Modelos Mistos (EMM)

umd <- t(Xc) %*% Rid %*% Xc
doisd<- t(Zc) %*% Rid %*% Xc
tresd<- t(Xc) %*% Rid %*% Zc
quatrod<- (t(Zc) %*% Rid %*% Zc) + Gid
cincod<- t(Xc) %*% Rid %*% Yc
seisd<- t(Zc) %*% Rid %*% Yc

LHSd <-cbind(rbind(umd,doisd),rbind(tresd, quatrod))
RHSd <-rbind(cincod, seisd)
paramd <-MASS::ginv(LHSd) %*% RHSd
paramd
##               [,1]
##  [1,] 2390.8381913
##  [2,]  213.9445733
##  [3,]  -21.3617856
##  [4,]   35.9815572
##  [5,]    4.0526622
##  [6,]    3.3563192
##  [7,]   30.3887522
##  [8,]  -20.7714030
##  [9,]    0.7236128
## [10,]  -10.0238951
## [11,]   -0.3445578
## [12,]    0.2125166
## [13,]    7.0579772
## [14,]   -7.0249981
## [15,]  -10.7482645
## [16,]   -8.8866313

Refinando…

paramd<- as.data.frame(paramd)

paramd<- tibble::rownames_to_column(paramd, var = "Parametros")
paramd$Parametros <- c("Média P42", "Média Ovos", " Sexo", "Covariável", "A1", "A1", "A2", "A2",
                    "A3", "A3", "A4", "A4", "A5", "A5", "A6", "A6")
colnames(paramd)[2] <- "Efeitos"
paramd
##    Parametros      Efeitos
## 1   Média P42 2390.8381913
## 2  Média Ovos  213.9445733
## 3        Sexo  -21.3617856
## 4  Covariável   35.9815572
## 5          A1    4.0526622
## 6          A1    3.3563192
## 7          A2   30.3887522
## 8          A2  -20.7714030
## 9          A3    0.7236128
## 10         A3  -10.0238951
## 11         A4   -0.3445578
## 12         A4    0.2125166
## 13         A5    7.0579772
## 14         A5   -7.0249981
## 15         A6  -10.7482645
## 16         A6   -8.8866313