Montar o vetor y; as matrizes X, Z, A, R e G e obter as soluções para as equações de modelos mistos do exemplo abaixo:
dados <- data.frame (Animal = c(1:12),
Pai = c ("A", "A", "A", "A", "A", "B", "B","B", "B","B", "B", "B"),
Mae = c ("C", "C", "D", "D", "D", "E", "E","E", "E","F", "F", "F"),
Incubacao = c(1,2,1,1,2,1,1,2,2,1,1,2),
Sexo = c(2,2,1,2,2,1,1,2,2,1,2,2),
PNascer = c(50,51,51,52,48,50,52,50,52,54,52,50),
P42 = c (2400,2300,2300,2400,2200,2250,2300,2200,2150,2350,2150,2300),
Ovos = c(182,180,0,181,182,0,0,178,175,0,179,183))dados## Animal Pai Mae Incubacao Sexo PNascer P42 Ovos
## 1 1 A C 1 2 50 2400 182
## 2 2 A C 2 2 51 2300 180
## 3 3 A D 1 1 51 2300 0
## 4 4 A D 1 2 52 2400 181
## 5 5 A D 2 2 48 2200 182
## 6 6 B E 1 1 50 2250 0
## 7 7 B E 1 1 52 2300 0
## 8 8 B E 2 2 50 2200 178
## 9 9 B E 2 2 52 2150 175
## 10 10 B F 1 1 54 2350 0
## 11 11 B F 1 2 52 2150 179
## 12 12 B F 2 2 50 2300 183Modelo 1(Característica 1= Peso aos 42 dias): Yijk = média + Incubacãoi + Sexoj + b *(PNascerijk - PNascerMédio) + aijk + eijk
P_42 <- as.matrix (dados$P42)
P_42## [,1]
## [1,] 2400
## [2,] 2300
## [3,] 2300
## [4,] 2400
## [5,] 2200
## [6,] 2250
## [7,] 2300
## [8,] 2200
## [9,] 2150
## [10,] 2350
## [11,] 2150
## [12,] 2300media_PN <- mean(dados$PNascer)
media_PN## [1] 51cov <- dados$PNascer - media_PN
cov## [1] -1 0 0 1 -3 -1 1 -1 1 3 1 -1X <- data.frame (media = c(rep(1,12)),
ind1= c( 1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, -1),
sex1= c(-1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1),
cov)
X <- as.matrix(X)
str(X)## num [1:12, 1:4] 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:4] "media" "ind1" "sex1" "cov"X## media ind1 sex1 cov
## [1,] 1 1 -1 -1
## [2,] 1 -1 -1 0
## [3,] 1 1 1 0
## [4,] 1 1 -1 1
## [5,] 1 -1 -1 -3
## [6,] 1 1 1 -1
## [7,] 1 1 1 1
## [8,] 1 -1 -1 -1
## [9,] 1 -1 -1 1
## [10,] 1 1 1 3
## [11,] 1 1 -1 1
## [12,] 1 -1 -1 -1Z <- as.matrix(diag(12))
Z## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 1 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 1 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 1 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 1P <- data.frame(Ind = c(21:26,1:12),
Par1 = c(0,0,0,0,0,0,21,21,21,21,21,22,22,22,22,22,22,22),
Par2 = c(0,0,0,0,0,0,23,23,24,24,24,25,25,25,25,26,26,26))
P## Ind Par1 Par2
## 1 21 0 0
## 2 22 0 0
## 3 23 0 0
## 4 24 0 0
## 5 25 0 0
## 6 26 0 0
## 7 1 21 23
## 8 2 21 23
## 9 3 21 24
## 10 4 21 24
## 11 5 21 24
## 12 6 22 25
## 13 7 22 25
## 14 8 22 25
## 15 9 22 25
## 16 10 22 26
## 17 11 22 26
## 18 12 22 26install.packages("AGHmatrix", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/ferre/OneDrive/Documentos/R/win-library/4.0'
## (as 'lib' is unspecified)## package 'AGHmatrix' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\ferre\AppData\Local\Temp\RtmpucgVpJ\downloaded_packageslibrary("AGHmatrix")## Warning: package 'AGHmatrix' was built under R version 4.0.5A <- 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 minutesA <-as.matrix (A[7:18,7:18])
A## 1 2 3 4 5 6 7 8 9 10 11 12
## 1 1.00 0.50 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 2 0.50 1.00 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 3 0.25 0.25 1.00 0.50 0.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 4 0.25 0.25 0.50 1.00 0.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 5 0.25 0.25 0.50 0.50 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 6 0.00 0.00 0.00 0.00 0.00 1.00 0.50 0.50 0.50 0.25 0.25 0.25
## 7 0.00 0.00 0.00 0.00 0.00 0.50 1.00 0.50 0.50 0.25 0.25 0.25
## 8 0.00 0.00 0.00 0.00 0.00 0.50 0.50 1.00 0.50 0.25 0.25 0.25
## 9 0.00 0.00 0.00 0.00 0.00 0.50 0.50 0.50 1.00 0.25 0.25 0.25
## 10 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 1.00 0.50 0.50
## 11 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.50 1.00 0.50
## 12 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.50 0.50 1.00I_P42 <- as.matrix (diag(12))
varE_P42 = 10000
R_P42 <- (I_P42* varE_P42)
R_P42## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 10000 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 10000 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 10000 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 10000 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 10000 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 10000 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 10000 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 10000 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 10000 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 10000 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 10000 0
## [12,] 0 0 0 0 0 0 0 0 0 0 0 10000Ri <- as.matrix (solve(R_P42))varG_P42 = 10000
G_P42 <- A * varG_P42
G_P42## 1 2 3 4 5 6 7 8 9 10 11 12
## 1 10000 5000 2500 2500 2500 0 0 0 0 0 0 0
## 2 5000 10000 2500 2500 2500 0 0 0 0 0 0 0
## 3 2500 2500 10000 5000 5000 0 0 0 0 0 0 0
## 4 2500 2500 5000 10000 5000 0 0 0 0 0 0 0
## 5 2500 2500 5000 5000 10000 0 0 0 0 0 0 0
## 6 0 0 0 0 0 10000 5000 5000 5000 2500 2500 2500
## 7 0 0 0 0 0 5000 10000 5000 5000 2500 2500 2500
## 8 0 0 0 0 0 5000 5000 10000 5000 2500 2500 2500
## 9 0 0 0 0 0 5000 5000 5000 10000 2500 2500 2500
## 10 0 0 0 0 0 2500 2500 2500 2500 10000 5000 5000
## 11 0 0 0 0 0 2500 2500 2500 2500 5000 10000 5000
## 12 0 0 0 0 0 2500 2500 2500 2500 5000 5000 10000Gi<- as.matrix(solve(G_P42))um<- t(X) %*% Ri %*% X
dois<- t(Z) %*% Ri %*% X
tres<- t(X) %*% Ri %*% Z
quatro<- (t(Z) %*% Ri %*% Z) + Gi
cinco<- t(X) %*% Ri %*% P_42
seis<- t(Z) %*% Ri %*% P_42
LHS <-cbind(rbind(um,dois),rbind(tres, quatro))
RHS <-rbind(cinco, seis)
install.packages("MASS", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/ferre/OneDrive/Documentos/R/win-library/4.0'
## (as 'lib' is unspecified)## package 'MASS' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'MASS'## Warning in file.copy(savedcopy, lib, recursive = TRUE): problem copying C:
## \Users\ferre\OneDrive\Documentos\R\win-library\4.0\00LOCK\MASS\libs\x64\MASS.dll
## to C:\Users\ferre\OneDrive\Documentos\R\win-library\4.0\MASS\libs\x64\MASS.dll:
## Permission denied## Warning: restored 'MASS'##
## The downloaded binary packages are in
## C:\Users\ferre\AppData\Local\Temp\RtmpucgVpJ\downloaded_packageslibrary(MASS)## Warning: package 'MASS' was built under R version 4.0.5param <-MASS::ginv(LHS) %*% RHS
install.packages("dplyr", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/ferre/OneDrive/Documentos/R/win-library/4.0'
## (as 'lib' is unspecified)## package 'dplyr' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'dplyr'## Warning in file.copy(savedcopy, lib, recursive = TRUE):
## problem copying C:\Users\ferre\OneDrive\Documentos\R\win-
## library\4.0\00LOCK\dplyr\libs\x64\dplyr.dll to C:
## \Users\ferre\OneDrive\Documentos\R\win-library\4.0\dplyr\libs\x64\dplyr.dll:
## Permission denied## Warning: restored 'dplyr'##
## The downloaded binary packages are in
## C:\Users\ferre\AppData\Local\Temp\RtmpucgVpJ\downloaded_packageslibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.0.5##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionparam2<- as.data.frame(param)
param2<- tibble::rownames_to_column(param2, var = "Efeitos")
param2$Efeitos <- c("Média", "Incubação", "Sexo", "P_Nascer", "A1", "A2", "A3", "A4", "A5", "A6", "A7", "A8", "A9", "A10", "A11", "A12")
colnames(param2)[2] <- "Peso 42 dias"
param2## Efeitos Peso 42 dias
## 1 Média 2276.013506
## 2 Incubação 32.035004
## 3 Sexo 3.375797
## 4 P_Nascer 5.677223
## 5 A1 55.284346
## 6 A2 41.415275
## 7 A3 7.573159
## 8 A4 41.264616
## 9 A5 3.524250
## 10 A6 -38.785451
## 11 A7 -25.903600
## 12 A8 -31.844917
## 13 A9 -52.296399
## 14 A10 -5.613008
## 15 A11 -66.244328
## 16 A12 8.897157Modelo 2(Característica 2 = Número de Ovos: Yijk= média + Ii + aij + Eij
N_ovos<- as.matrix (dados$Ovos)
N_ovos## [,1]
## [1,] 182
## [2,] 180
## [3,] 0
## [4,] 181
## [5,] 182
## [6,] 0
## [7,] 0
## [8,] 178
## [9,] 175
## [10,] 0
## [11,] 179
## [12,] 183incb <- dados$Incubacao
incb [incb == 2 ] <- -1
Xb <- data.frame (media = c(rep (1,12)),
Incb1= c(incb))
Xb <- as.matrix(Xb)
Xb## media Incb1
## [1,] 1 1
## [2,] 1 -1
## [3,] 1 1
## [4,] 1 1
## [5,] 1 -1
## [6,] 1 1
## [7,] 1 1
## [8,] 1 -1
## [9,] 1 -1
## [10,] 1 1
## [11,] 1 1
## [12,] 1 -1Zb <- as.matrix(diag(12))
Zb## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 1 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 1 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 1 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 1A## 1 2 3 4 5 6 7 8 9 10 11 12
## 1 1.00 0.50 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 2 0.50 1.00 0.25 0.25 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 3 0.25 0.25 1.00 0.50 0.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 4 0.25 0.25 0.50 1.00 0.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 5 0.25 0.25 0.50 0.50 1.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## 6 0.00 0.00 0.00 0.00 0.00 1.00 0.50 0.50 0.50 0.25 0.25 0.25
## 7 0.00 0.00 0.00 0.00 0.00 0.50 1.00 0.50 0.50 0.25 0.25 0.25
## 8 0.00 0.00 0.00 0.00 0.00 0.50 0.50 1.00 0.50 0.25 0.25 0.25
## 9 0.00 0.00 0.00 0.00 0.00 0.50 0.50 0.50 1.00 0.25 0.25 0.25
## 10 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 1.00 0.50 0.50
## 11 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.50 1.00 0.50
## 12 0.00 0.00 0.00 0.00 0.00 0.25 0.25 0.25 0.25 0.50 0.50 1.00I_ovos <- as.matrix (diag(12))
varE_ovos = 180
R_ovos <- (I_ovos* varE_ovos)
R_ovos## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 180 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 180 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 180 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 180 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 180 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 180 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 180 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 180 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 180 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 180 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 180 0
## [12,] 0 0 0 0 0 0 0 0 0 0 0 180Rbi <- as.matrix (solve(R_ovos))
varG_ovos = 50
G_ovos <- A * varG_ovos
G_ovos## 1 2 3 4 5 6 7 8 9 10 11 12
## 1 50.0 25.0 12.5 12.5 12.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## 2 25.0 50.0 12.5 12.5 12.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## 3 12.5 12.5 50.0 25.0 25.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## 4 12.5 12.5 25.0 50.0 25.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## 5 12.5 12.5 25.0 25.0 50.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## 6 0.0 0.0 0.0 0.0 0.0 50.0 25.0 25.0 25.0 12.5 12.5 12.5
## 7 0.0 0.0 0.0 0.0 0.0 25.0 50.0 25.0 25.0 12.5 12.5 12.5
## 8 0.0 0.0 0.0 0.0 0.0 25.0 25.0 50.0 25.0 12.5 12.5 12.5
## 9 0.0 0.0 0.0 0.0 0.0 25.0 25.0 25.0 50.0 12.5 12.5 12.5
## 10 0.0 0.0 0.0 0.0 0.0 12.5 12.5 12.5 12.5 50.0 25.0 25.0
## 11 0.0 0.0 0.0 0.0 0.0 12.5 12.5 12.5 12.5 25.0 50.0 25.0
## 12 0.0 0.0 0.0 0.0 0.0 12.5 12.5 12.5 12.5 25.0 25.0 50.0Gbi<- as.matrix(solve(G_ovos))Resumindo LHS = parâmetros * RHS . Logo, parâmetros = inversa de LHS * RHS
um<- t(Xb) %*% Rbi %*% Xb
dois<- t(Zb) %*% Rbi %*% Xb
tres<- t(Xb) %*% Rbi %*% Zb
quatro<- (t(Zb) %*% Rbi %*% Zb) + Gbi
cinco<- t(Xb) %*% Rbi %*% N_ovos
seis<- t(Z) %*% Ri %*% N_ovos
LHS <-cbind(rbind(um,dois),rbind(tres, quatro))
RHS <-rbind(cinco, seis)
param <-MASS::ginv(LHS) %*% RHS
param2<- as.data.frame(param)
param2<- tibble::rownames_to_column(param2, var = "Efeitos")
param2$Efeitos <- c("Média", "Incubação", "A1", "A2", "A3", "A4", "A5", "A6", "A7", "A8", "A9", "A10", "A11", "A12")
colnames(param2)[2] <- "Número de Ovos"
param2## Efeitos Número de Ovos
## 1 Média 221.86550
## 2 Incubação -58.61401
## 3 A1 -74.26642
## 4 A2 -88.56691
## 5 A3 -79.23139
## 6 A4 -78.83407
## 7 A5 -93.12798
## 8 A6 -96.58024
## 9 A7 -96.58024
## 10 A8 -110.48561
## 11 A9 -110.49220
## 12 A10 -87.83048
## 13 A11 -87.43755
## 14 A12 -101.72487Modelo: Yijk = Sexoi + Aij + Eij
Passos: i) Instalar o pacotesreshape2 ii) 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”: iii) Ordenar o dados_melt por animal
dados_melt <- reshape2::melt (dados, id.vars=c("Animal", "Pai", "Mae", "Incubacao", "Sexo", "PNascer"))
dadosc<- dados_melt[order(dados_melt$Animal),]
dadosc## Animal Pai Mae Incubacao Sexo PNascer variable value
## 1 1 A C 1 2 50 P42 2400
## 13 1 A C 1 2 50 Ovos 182
## 2 2 A C 2 2 51 P42 2300
## 14 2 A C 2 2 51 Ovos 180
## 3 3 A D 1 1 51 P42 2300
## 15 3 A D 1 1 51 Ovos 0
## 4 4 A D 1 2 52 P42 2400
## 16 4 A D 1 2 52 Ovos 181
## 5 5 A D 2 2 48 P42 2200
## 17 5 A D 2 2 48 Ovos 182
## 6 6 B E 1 1 50 P42 2250
## 18 6 B E 1 1 50 Ovos 0
## 7 7 B E 1 1 52 P42 2300
## 19 7 B E 1 1 52 Ovos 0
## 8 8 B E 2 2 50 P42 2200
## 20 8 B E 2 2 50 Ovos 178
## 9 9 B E 2 2 52 P42 2150
## 21 9 B E 2 2 52 Ovos 175
## 10 10 B F 1 1 54 P42 2350
## 22 10 B F 1 1 54 Ovos 0
## 11 11 B F 1 2 52 P42 2150
## 23 11 B F 1 2 52 Ovos 179
## 12 12 B F 2 2 50 P42 2300
## 24 12 B F 2 2 50 Ovos 183##Passo 2: Escrever o vetor Y
Yc <-as.matrix(dadosc$value)
Yc## [,1]
## [1,] 2400
## [2,] 182
## [3,] 2300
## [4,] 180
## [5,] 2300
## [6,] 0
## [7,] 2400
## [8,] 181
## [9,] 2200
## [10,] 182
## [11,] 2250
## [12,] 0
## [13,] 2300
## [14,] 0
## [15,] 2200
## [16,] 178
## [17,] 2150
## [18,] 175
## [19,] 2350
## [20,] 0
## [21,] 2150
## [22,] 179
## [23,] 2300
## [24,] 183a1 <- dadosc$variable =="P42" & dadosc$Sexo == 1
a1 <- as.numeric(a1)
a2 <- dadosc$variable =="Ovos" & dadosc$Sexo == 1
a2 <- as.numeric(a2)
a3 <- dadosc$variable =="P42" & dadosc$Sexo == 2
a3 <- as.numeric(a3)
a4 <- dadosc$variable =="Ovos" & dadosc$Sexo == 2
a4 <- as.numeric(a4)
Xc <- cbind(a1, a2, a3, a4)
Xc## a1 a2 a3 a4
## [1,] 0 0 1 0
## [2,] 0 0 0 1
## [3,] 0 0 1 0
## [4,] 0 0 0 1
## [5,] 1 0 0 0
## [6,] 0 1 0 0
## [7,] 0 0 1 0
## [8,] 0 0 0 1
## [9,] 0 0 1 0
## [10,] 0 0 0 1
## [11,] 1 0 0 0
## [12,] 0 1 0 0
## [13,] 1 0 0 0
## [14,] 0 1 0 0
## [15,] 0 0 1 0
## [16,] 0 0 0 1
## [17,] 0 0 1 0
## [18,] 0 0 0 1
## [19,] 1 0 0 0
## [20,] 0 1 0 0
## [21,] 0 0 1 0
## [22,] 0 0 0 1
## [23,] 0 0 1 0
## [24,] 0 0 0 1Zc <- as.matrix(diag(24))
Zc## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 1 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 1 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 1 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 1 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 1 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 1 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 1 0 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 1 0 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 1 0 0
## [12,] 0 0 0 0 0 0 0 0 0 0 0 1 0
## [13,] 0 0 0 0 0 0 0 0 0 0 0 0 1
## [14,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [15,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [16,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [17,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [18,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [19,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [20,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [21,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [22,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [23,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [24,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 0
## [12,] 0 0 0 0 0 0 0 0 0 0 0
## [13,] 0 0 0 0 0 0 0 0 0 0 0
## [14,] 1 0 0 0 0 0 0 0 0 0 0
## [15,] 0 1 0 0 0 0 0 0 0 0 0
## [16,] 0 0 1 0 0 0 0 0 0 0 0
## [17,] 0 0 0 1 0 0 0 0 0 0 0
## [18,] 0 0 0 0 1 0 0 0 0 0 0
## [19,] 0 0 0 0 0 1 0 0 0 0 0
## [20,] 0 0 0 0 0 0 1 0 0 0 0
## [21,] 0 0 0 0 0 0 0 1 0 0 0
## [22,] 0 0 0 0 0 0 0 0 1 0 0
## [23,] 0 0 0 0 0 0 0 0 0 1 0
## [24,] 0 0 0 0 0 0 0 0 0 0 1Ac <- A
Ac <- as.matrix(Ac)Ic <- as.matrix (diag(12))
Rzero <- rbind(c (10000, -250), c (-250, 180))
Rc <- Ic %x% Rzero
Rc## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 10000 -250 0 0 0 0 0 0 0 0 0 0 0
## [2,] -250 180 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 10000 -250 0 0 0 0 0 0 0 0 0
## [4,] 0 0 -250 180 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 10000 -250 0 0 0 0 0 0 0
## [6,] 0 0 0 0 -250 180 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 10000 -250 0 0 0 0 0
## [8,] 0 0 0 0 0 0 -250 180 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 10000 -250 0 0 0
## [10,] 0 0 0 0 0 0 0 0 -250 180 0 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 10000 -250 0
## [12,] 0 0 0 0 0 0 0 0 0 0 -250 180 0
## [13,] 0 0 0 0 0 0 0 0 0 0 0 0 10000
## [14,] 0 0 0 0 0 0 0 0 0 0 0 0 -250
## [15,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [16,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [17,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [18,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [19,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [20,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [21,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [22,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [23,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [24,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 0
## [12,] 0 0 0 0 0 0 0 0 0 0 0
## [13,] -250 0 0 0 0 0 0 0 0 0 0
## [14,] 180 0 0 0 0 0 0 0 0 0 0
## [15,] 0 10000 -250 0 0 0 0 0 0 0 0
## [16,] 0 -250 180 0 0 0 0 0 0 0 0
## [17,] 0 0 0 10000 -250 0 0 0 0 0 0
## [18,] 0 0 0 -250 180 0 0 0 0 0 0
## [19,] 0 0 0 0 0 10000 -250 0 0 0 0
## [20,] 0 0 0 0 0 -250 180 0 0 0 0
## [21,] 0 0 0 0 0 0 0 10000 -250 0 0
## [22,] 0 0 0 0 0 0 0 -250 180 0 0
## [23,] 0 0 0 0 0 0 0 0 0 10000 -250
## [24,] 0 0 0 0 0 0 0 0 0 -250 180Ric <- solve(Rc)
Ric <- as.matrix(Ric)
Gzero <- rbind(c (10000, -100), c (-100, 50))
Gzero <-as.matrix(Gzero)
Gc <- Ac %x% Gzero
Gc## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## [1,] 10000 -100.0 5000 -50.0 2500 -25.0 2500 -25.0 2500 -25.0 0
## [2,] -100 50.0 -50 25.0 -25 12.5 -25 12.5 -25 12.5 0
## [3,] 5000 -50.0 10000 -100.0 2500 -25.0 2500 -25.0 2500 -25.0 0
## [4,] -50 25.0 -100 50.0 -25 12.5 -25 12.5 -25 12.5 0
## [5,] 2500 -25.0 2500 -25.0 10000 -100.0 5000 -50.0 5000 -50.0 0
## [6,] -25 12.5 -25 12.5 -100 50.0 -50 25.0 -50 25.0 0
## [7,] 2500 -25.0 2500 -25.0 5000 -50.0 10000 -100.0 5000 -50.0 0
## [8,] -25 12.5 -25 12.5 -50 25.0 -100 50.0 -50 25.0 0
## [9,] 2500 -25.0 2500 -25.0 5000 -50.0 5000 -50.0 10000 -100.0 0
## [10,] -25 12.5 -25 12.5 -50 25.0 -50 25.0 -100 50.0 0
## [11,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 10000
## [12,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 -100
## [13,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 5000
## [14,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 -50
## [15,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 5000
## [16,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 -50
## [17,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 5000
## [18,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 -50
## [19,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 2500
## [20,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 -25
## [21,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 2500
## [22,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 -25
## [23,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 2500
## [24,] 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0 -25
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## [1,] 0.0 0 0.0 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 0.0 0 0.0
## [3,] 0.0 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0
## [4,] 0.0 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0
## [5,] 0.0 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0
## [6,] 0.0 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0
## [7,] 0.0 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0
## [8,] 0.0 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0
## [9,] 0.0 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0
## [10,] 0.0 0 0.0 0 0.0 0 0.0 0 0.0 0 0.0
## [11,] -100.0 5000 -50.0 5000 -50.0 5000 -50.0 2500 -25.0 2500 -25.0
## [12,] 50.0 -50 25.0 -50 25.0 -50 25.0 -25 12.5 -25 12.5
## [13,] -50.0 10000 -100.0 5000 -50.0 5000 -50.0 2500 -25.0 2500 -25.0
## [14,] 25.0 -100 50.0 -50 25.0 -50 25.0 -25 12.5 -25 12.5
## [15,] -50.0 5000 -50.0 10000 -100.0 5000 -50.0 2500 -25.0 2500 -25.0
## [16,] 25.0 -50 25.0 -100 50.0 -50 25.0 -25 12.5 -25 12.5
## [17,] -50.0 5000 -50.0 5000 -50.0 10000 -100.0 2500 -25.0 2500 -25.0
## [18,] 25.0 -50 25.0 -50 25.0 -100 50.0 -25 12.5 -25 12.5
## [19,] -25.0 2500 -25.0 2500 -25.0 2500 -25.0 10000 -100.0 5000 -50.0
## [20,] 12.5 -25 12.5 -25 12.5 -25 12.5 -100 50.0 -50 25.0
## [21,] -25.0 2500 -25.0 2500 -25.0 2500 -25.0 5000 -50.0 10000 -100.0
## [22,] 12.5 -25 12.5 -25 12.5 -25 12.5 -50 25.0 -100 50.0
## [23,] -25.0 2500 -25.0 2500 -25.0 2500 -25.0 5000 -50.0 5000 -50.0
## [24,] 12.5 -25 12.5 -25 12.5 -25 12.5 -50 25.0 -50 25.0
## [,23] [,24]
## [1,] 0 0.0
## [2,] 0 0.0
## [3,] 0 0.0
## [4,] 0 0.0
## [5,] 0 0.0
## [6,] 0 0.0
## [7,] 0 0.0
## [8,] 0 0.0
## [9,] 0 0.0
## [10,] 0 0.0
## [11,] 2500 -25.0
## [12,] -25 12.5
## [13,] 2500 -25.0
## [14,] -25 12.5
## [15,] 2500 -25.0
## [16,] -25 12.5
## [17,] 2500 -25.0
## [18,] -25 12.5
## [19,] 5000 -50.0
## [20,] -50 25.0
## [21,] 5000 -50.0
## [22,] -50 25.0
## [23,] 10000 -100.0
## [24,] -100 50.0Gic <- solve(Gc)
Gic <- as.matrix(Gic)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,] 2318.35681119
## [2,] 0.16552665
## [3,] 2262.50000000
## [4,] 180.00000000
## [5,] 70.66193506
## [6,] 0.33515386
## [7,] 36.32357660
## [8,] 0.18607513
## [9,] 5.94025487
## [10,] 0.26466787
## [11,] 58.96857122
## [12,] 0.26164722
## [13,] -8.45186429
## [14,] 0.56650485
## [15,] -47.95171744
## [16,] -0.38626594
## [17,] -31.03379450
## [18,] -0.43232959
## [19,] -46.43093879
## [20,] -0.61290468
## [21,] -64.10263058
## [22,] -0.92865007
## [23,] -0.38198768
## [24,] -0.10817894
## [25,] -49.36372159
## [26,] -0.07602371
## [27,] 2.39507238
## [28,] 0.26819740Refinando…
paramc<- as.data.frame(paramc)
paramc<- tibble::rownames_to_column(paramc, var = "Parametros")
paramc$Parametros <- c("M_P42", "M_Ovos", "F_P42", "F_Ovos", "A1", "A1", "A2", "A2",
"A3", "A3", "A4", "A4", "A5", "A5", "A6", "A6", "A7", "A7",
"A8", "A8", "A9", "A9", "A10", "A10", "A11", "A11", "A12", "A12")
colnames(paramc)[2] <- "Efeitos"
paramc## Parametros Efeitos
## 1 M_P42 2318.35681119
## 2 M_Ovos 0.16552665
## 3 F_P42 2262.50000000
## 4 F_Ovos 180.00000000
## 5 A1 70.66193506
## 6 A1 0.33515386
## 7 A2 36.32357660
## 8 A2 0.18607513
## 9 A3 5.94025487
## 10 A3 0.26466787
## 11 A4 58.96857122
## 12 A4 0.26164722
## 13 A5 -8.45186429
## 14 A5 0.56650485
## 15 A6 -47.95171744
## 16 A6 -0.38626594
## 17 A7 -31.03379450
## 18 A7 -0.43232959
## 19 A8 -46.43093879
## 20 A8 -0.61290468
## 21 A9 -64.10263058
## 22 A9 -0.92865007
## 23 A10 -0.38198768
## 24 A10 -0.10817894
## 25 A11 -49.36372159
## 26 A11 -0.07602371
## 27 A12 2.39507238
## 28 A12 0.26819740