Los datos de las mujeres fueron tomados de Libro: Epidemiologia y estadística para principiantes pagina 216 y los de los hombre fueron creados - pseudo reales.
Supongamos que tenemos una nueva prueba para detectar alfa-fetoproteína como predictor de un defecto del cierre del tubo neural que posee una sensibilidad del 90% y una especificidad del 95%.
Sexo: Mujer - hombre -> 2 niveles, riesgo: alto y bajo -> 2 niveles
alfa-fetoproteína: Anormal, normal Defecto del tubo neural Sano
| EVENTOS AL NACER | |||||||
|---|---|---|---|---|---|---|---|
| Poblacion | Sexo | Riesgo | alfa-fetoproteína | Defecto del tubo neuronal | Sano | Total | |
| 1 | Mujer | Alto | Anormal | 87 | 18 | ||
| Normal | 13 | 9882 | 10 000 | ||||
| 2 | Mujer | Bajo | Anormal | 128 | 179 | ||
| Normal | 19 | 99 674 | 100 000 | ||||
| 3 | Hombre | Alto | Anormal | 75 | 116 | ||
| Normal | 86 | 9723 | 10 000 | ||||
| 4 | Hombre | Bajo | Anormal | 108 | 267 | ||
| Normal | 61 | 99 564 | 100 000 |
Reescribimos la tabla anterior
pob1 <- c(9882,18,13,87)
pob2 <- c(99674,179,19,128)
pob3 <- c(9723,116,86,75)
pob4 <- c(99564,267,61,108)
Total <- c(sum(pob1), sum(pob2), sum(pob3), sum(pob4))
pob <- data.frame(t(data.frame(pob1, pob2, pob3, pob4)), Total,
row.names = c("Poblacion 1", "Poblacion 2","Poblacion 3","Poblacion 4"))
library(knitr)
knitr::kable(pob, col.names = c("C+P+", "C+P-", "C-P+", "C-P-", "Total"),
row.names = T)pob1 <- c(87,13,18,9882)
pob2 <- c(128,19,179,99674)
pob3 <- c(75,86,116,9723)
pob4 <- c(108,61,267,99564)
Total <- c(sum(pob1), sum(pob2), sum(pob3), sum(pob4))
pob <- data.frame(t(data.frame(pob1, pob2, pob3, pob4)), Total,
row.names = c("Poblacion 1", "Poblacion 2","Poblacion 3","Poblacion 4"))
library(knitr)
knitr::kable(pob, col.names = c("C+P+", "C+P-", "C-P+", "C-P-", "Total"),
row.names = T)| C+P+ | C+P- | C-P+ | C-P- | Total | |
|---|---|---|---|---|---|
| Poblacion 1 | 87 | 13 | 18 | 9882 | 1e+04 |
| Poblacion 2 | 128 | 19 | 179 | 99674 | 1e+05 |
| Poblacion 3 | 75 | 86 | 116 | 9723 | 1e+04 |
| Poblacion 4 | 108 | 61 | 267 | 99564 | 1e+05 |
La distribución de la tabla anterior viene dada por:
pi_pob1 <- pob1/sum(pob1)
pi_pob2 <- pob2/sum(pob2)
pi_pob3 <- pob3/sum(pob3)
pi_pob4 <- pob1/sum(pob4)
Total_prob <- c(1,1,1,1)
prob_pob <- data.frame(t(data.frame(pi_pob1, pi_pob2, pi_pob3, pi_pob4)), Total_prob,
row.names = c("Poblacion 1", "Poblacion 2","Poblacion 3","Poblacion 4"))
library(knitr)
knitr::kable(prob_pob, col.names = c("C+P+", "C+P-", "C-P+", "C-P-", "Total"),
row.names = T)| C+P+ | C+P- | C-P+ | C-P- | Total | |
|---|---|---|---|---|---|
| Poblacion 1 | 0.00870 | 0.00130 | 0.00180 | 0.98820 | 1 |
| Poblacion 2 | 0.00128 | 0.00019 | 0.00179 | 0.99674 | 1 |
| Poblacion 3 | 0.00750 | 0.00860 | 0.01160 | 0.97230 | 1 |
| Poblacion 4 | 0.00087 | 0.00013 | 0.00018 | 0.09882 | 1 |
pob <- c(pob1,pob2,pob3,pob4)
N <- pob
pi_pob1 <- pob1/sum(pob1)
pi_pob2 <- pob2/sum(pob2)
pi_pob3 <- pob3/sum(pob3)
pi_pob4 <- pob1/sum(pob4)
pi <- c(pi_pob1, pi_pob2, pi_pob3, pi_pob4)A <- matrix(c(1, 0, 0, 0,
1, 1, 0, 0,
0, 0, 0, 1,
0, 0, 1, 1),4, byrow = T)
identidad <- function(k)diag(rep(1,k))
Nro.Estratos <- 4
A <- kronecker(identidad(Nro.Estratos),A)
A## [,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,] 1 1 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 1 0 0 0 0 0 0 0 0 0
## [4,] 0 0 1 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 1 1 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 1 0 0 0 0 0
## [8,] 0 0 0 0 0 0 1 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 1 1 0 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 0 1 0
## [12,] 0 0 0 0 0 0 0 0 0 0 1 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 1
## [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
## [,14] [,15] [,16]
## [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,] 0 0 0
## [12,] 0 0 0
## [13,] 0 0 0
## [14,] 1 0 0
## [15,] 0 0 1
## [16,] 0 1 1A_Pi <- A%*%pi
A_Pi## [,1]
## [1,] 0.00870
## [2,] 0.01000
## [3,] 0.98820
## [4,] 0.99000
## [5,] 0.00128
## [6,] 0.00147
## [7,] 0.99674
## [8,] 0.99853
## [9,] 0.00750
## [10,] 0.01610
## [11,] 0.97230
## [12,] 0.98390
## [13,] 0.00087
## [14,] 0.00100
## [15,] 0.09882
## [16,] 0.09900LogNatu <- log(A_Pi)
LogNatu## [,1]
## [1,] -4.744432253
## [2,] -4.605170186
## [3,] -0.011870173
## [4,] -0.010050336
## [5,] -6.660895201
## [6,] -6.522492878
## [7,] -0.003265325
## [8,] -0.001471082
## [9,] -4.892852258
## [10,] -4.128936007
## [11,] -0.028090880
## [12,] -0.016231013
## [13,] -7.047017346
## [14,] -6.907755279
## [15,] -2.314455266
## [16,] -2.312635429k <- matrix(c(1, -1, 0, 0,
0, 0, 1, -1),2,4, byrow = T)
k <- kronecker(identidad(Nro.Estratos),k)
k## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1 -1 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 1 -1 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 1 -1 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 1 -1 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 1 -1 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 1 -1 0 0
## [7,] 0 0 0 0 0 0 0 0 0 0 0 0 1 -1
## [8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,15] [,16]
## [1,] 0 0
## [2,] 0 0
## [3,] 0 0
## [4,] 0 0
## [5,] 0 0
## [6,] 0 0
## [7,] 0 0
## [8,] 1 -1kln <- k%*%LogNatu
kln## [,1]
## [1,] -0.139262067
## [2,] -0.001819837
## [3,] -0.138402323
## [4,] -0.001794244
## [5,] -0.763916251
## [6,] -0.011859867
## [7,] -0.139262067
## [8,] -0.001819837ekln <- exp(kln)
ekln## [,1]
## [1,] 0.8700000
## [2,] 0.9981818
## [3,] 0.8707483
## [4,] 0.9982074
## [5,] 0.4658385
## [6,] 0.9882102
## [7,] 0.8700000
## [8,] 0.9981818Q <- diag(1,2*Nro.Estratos)
Q## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0 0 0
## [3,] 0 0 1 0 0 0 0 0
## [4,] 0 0 0 1 0 0 0 0
## [5,] 0 0 0 0 1 0 0 0
## [6,] 0 0 0 0 0 1 0 0
## [7,] 0 0 0 0 0 0 1 0
## [8,] 0 0 0 0 0 0 0 1#Calcular f_exp-estimado
f <- Q%*%ekln
f## [,1]
## [1,] 0.8700000
## [2,] 0.9981818
## [3,] 0.8707483
## [4,] 0.9982074
## [5,] 0.4658385
## [6,] 0.9882102
## [7,] 0.8700000
## [8,] 0.9981818Para el modelo \(\hat { f } =X\beta+\epsilon\), donde \(Var\left( \epsilon\right)={\Sigma}_{\hat{\epsilon}}\) y \(\epsilon\sim AN \left( 0,{\Sigma}_{\hat{\epsilon}} \right)\) {AN asintóticamente normal} el estimador via mínimos cuadrados ponderados de \(\beta\) el cual es un estimador BAN (best asymptotic normal) es:
\[\begin{equation} \hat { \beta }={ \left( X'{ \hat { \Sigma } }_{ \hat { f } }^{ -1 }X \right) }^{ -1 }X'{ \hat { \Sigma } }_{ \hat { f } }^{ -1 }\hat { f } \end{equation}\]
Para estimar \(\beta\), comenzamos a hallar o definir:
########################## Matriz diseño
repita<-function(X,X2,k){
temp<-NULL
for(i in 1:nrow(X))temp<-rbind(temp,cbind
(matrix(rep(X[i,],k),ncol=ncol(X),byrow=T),X2))
return(temp)
}
nro.vars<-2
niveles<-c(2,2)
bloque <- function(k){rbind(0,identidad(k))}
matriz.diseno <- function(nro.vars,niveles){
niveles <- niveles-1
temp1 <- bloque(niveles[nro.vars])
for(i in (nro.vars-1):1){
temp2<-bloque(niveles[i])
temp1<-repita(temp2,temp1,nrow(temp1))
}
temp <- cbind(1,temp1)
temp <- temp[rep(1:nrow(temp),each=2),]
indi <- matrix(rep(c(0,1),nrow(temp)/2),ncol=1)
temp2 <- NULL
for(ii in 1:ncol(temp)) temp2<-cbind(temp2,indi)
return(cbind(temp,temp*temp2))
}
X <- matriz.diseno(nro.vars, niveles)
X## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 0 0
## [2,] 1 0 0 1 0 0
## [3,] 1 0 1 0 0 0
## [4,] 1 0 1 1 0 1
## [5,] 1 1 0 0 0 0
## [6,] 1 1 0 1 1 0
## [7,] 1 1 1 0 0 0
## [8,] 1 1 1 1 1 1N <- pob
crea.bloque <- function(mat1,mat2){
nf1<-nrow(mat1);nc1<-ncol(mat1)
nf2<-nrow(mat2);nc2<-ncol(mat2)
mat1<-cbind(mat1,matrix(0,nrow=nf1,ncol=nc2))
mat2<-cbind(matrix(0,nrow=nf2,ncol=nc1),mat2)
matdef<-rbind(mat1,mat2)
return(matdef)
}
donde.voy<-0
for(i in 1:Nro.Estratos){
print("====================")
identifica<-c("Estrato ",i)
print(identifica)
temp<-N[(donde.voy+1):(donde.voy+4)]
print("Frecuencias de la subpoblacion")
print(temp)
donde.voy<-donde.voy+4
n.temp<-sum(temp)
probab<-matrix(temp/n.temp,ncol=1)
temp1<- probab%*%t(probab)
temp2<-diag(as.vector(probab))
varcov.p<-(temp2-temp1)/n.temp
print("p estimado del estrato")
print(t(probab))
print("Matriz de varianzas y covarianzas estimada del estrato")
print(varcov.p)
if(i==1){
varcov.grande<-varcov.p
prob.grande<-probab
}
else{
varcov.grande <- crea.bloque(varcov.grande,varcov.p)
prob.grande <- rbind(prob.grande,probab) # Todas las probabilidades organizadas en columna
}
}## [1] "===================="
## [1] "Estrato " "1"
## [1] "Frecuencias de la subpoblacion"
## [1] 87 13 18 9882
## [1] "p estimado del estrato"
## [,1] [,2] [,3] [,4]
## [1,] 0.0087 0.0013 0.0018 0.9882
## [1] "Matriz de varianzas y covarianzas estimada del estrato"
## [,1] [,2] [,3] [,4]
## [1,] 8.62431e-07 -1.13100e-09 -1.56600e-09 -8.597340e-07
## [2,] -1.13100e-09 1.29831e-07 -2.34000e-10 -1.284660e-07
## [3,] -1.56600e-09 -2.34000e-10 1.79676e-07 -1.778760e-07
## [4,] -8.59734e-07 -1.28466e-07 -1.77876e-07 1.166076e-06
## [1] "===================="
## [1] "Estrato " "2"
## [1] "Frecuencias de la subpoblacion"
## [1] 128 19 179 99674
## [1] "p estimado del estrato"
## [,1] [,2] [,3] [,4]
## [1,] 0.00128 0.00019 0.00179 0.99674
## [1] "Matriz de varianzas y covarianzas estimada del estrato"
## [,1] [,2] [,3] [,4]
## [1,] 1.278362e-08 -2.432000e-12 -2.291200e-11 -1.275827e-08
## [2,] -2.432000e-12 1.899639e-09 -3.401000e-12 -1.893806e-09
## [3,] -2.291200e-11 -3.401000e-12 1.786796e-08 -1.784165e-08
## [4,] -1.275827e-08 -1.893806e-09 -1.784165e-08 3.249372e-08
## [1] "===================="
## [1] "Estrato " "3"
## [1] "Frecuencias de la subpoblacion"
## [1] 75 86 116 9723
## [1] "p estimado del estrato"
## [,1] [,2] [,3] [,4]
## [1,] 0.0075 0.0086 0.0116 0.9723
## [1] "Matriz de varianzas y covarianzas estimada del estrato"
## [,1] [,2] [,3] [,4]
## [1,] 7.44375e-07 -6.45000e-09 -8.700000e-09 -7.292250e-07
## [2,] -6.45000e-09 8.52604e-07 -9.976000e-09 -8.361780e-07
## [3,] -8.70000e-09 -9.97600e-09 1.146544e-06 -1.127868e-06
## [4,] -7.29225e-07 -8.36178e-07 -1.127868e-06 2.693271e-06
## [1] "===================="
## [1] "Estrato " "4"
## [1] "Frecuencias de la subpoblacion"
## [1] 108 61 267 99564
## [1] "p estimado del estrato"
## [,1] [,2] [,3] [,4]
## [1,] 0.00108 0.00061 0.00267 0.99564
## [1] "Matriz de varianzas y covarianzas estimada del estrato"
## [,1] [,2] [,3] [,4]
## [1,] 1.078834e-08 -6.588000e-12 -2.883600e-11 -1.075291e-08
## [2,] -6.588000e-12 6.096279e-09 -1.628700e-11 -6.073404e-09
## [3,] -2.883600e-11 -1.628700e-11 2.662871e-08 -2.658359e-08
## [4,] -1.075291e-08 -6.073404e-09 -2.658359e-08 4.340990e-08# Calcular (A*sigma pi-gorro*A")
var.cov.lineal <- A%*%varcov.grande%*%t(A)
var.cov.lineal## [,1] [,2] [,3] [,4] [,5]
## [1,] 8.62431e-07 8.613e-07 -8.597340e-07 -8.613e-07 0.000000e+00
## [2,] 8.61300e-07 9.900e-07 -9.882000e-07 -9.900e-07 0.000000e+00
## [3,] -8.59734e-07 -9.882e-07 1.166076e-06 9.882e-07 0.000000e+00
## [4,] -8.61300e-07 -9.900e-07 9.882000e-07 9.900e-07 0.000000e+00
## [5,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 1.278362e-08
## [6,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 1.278118e-08
## [7,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 -1.275827e-08
## [8,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 -1.278118e-08
## [9,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 0.000000e+00
## [10,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 0.000000e+00
## [11,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 0.000000e+00
## [12,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 0.000000e+00
## [13,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 0.000000e+00
## [14,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 0.000000e+00
## [15,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 0.000000e+00
## [16,] 0.00000e+00 0.000e+00 0.000000e+00 0.000e+00 0.000000e+00
## [,6] [,7] [,8] [,9] [,10]
## [1,] 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.000000e+00
## [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.000000e+00
## [3,] 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.000000e+00
## [4,] 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.000000e+00
## [5,] 1.278118e-08 -1.275827e-08 -1.278118e-08 0.00000e+00 0.000000e+00
## [6,] 1.467839e-08 -1.465208e-08 -1.467839e-08 0.00000e+00 0.000000e+00
## [7,] -1.465208e-08 3.249372e-08 1.465208e-08 0.00000e+00 0.000000e+00
## [8,] -1.467839e-08 1.465208e-08 1.467839e-08 0.00000e+00 0.000000e+00
## [9,] 0.000000e+00 0.000000e+00 0.000000e+00 7.44375e-07 7.379250e-07
## [10,] 0.000000e+00 0.000000e+00 0.000000e+00 7.37925e-07 1.584079e-06
## [11,] 0.000000e+00 0.000000e+00 0.000000e+00 -7.29225e-07 -1.565403e-06
## [12,] 0.000000e+00 0.000000e+00 0.000000e+00 -7.37925e-07 -1.584079e-06
## [13,] 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.000000e+00
## [14,] 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.000000e+00
## [15,] 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.000000e+00
## [16,] 0.000000e+00 0.000000e+00 0.000000e+00 0.00000e+00 0.000000e+00
## [,11] [,12] [,13] [,14] [,15]
## [1,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [3,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [4,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [5,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [6,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [7,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [8,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [9,] -7.292250e-07 -7.379250e-07 0.000000e+00 0.000000e+00 0.000000e+00
## [10,] -1.565403e-06 -1.584079e-06 0.000000e+00 0.000000e+00 0.000000e+00
## [11,] 2.693271e-06 1.565403e-06 0.000000e+00 0.000000e+00 0.000000e+00
## [12,] 1.565403e-06 1.584079e-06 0.000000e+00 0.000000e+00 0.000000e+00
## [13,] 0.000000e+00 0.000000e+00 1.078834e-08 1.078175e-08 -1.075291e-08
## [14,] 0.000000e+00 0.000000e+00 1.078175e-08 1.687144e-08 -1.682632e-08
## [15,] 0.000000e+00 0.000000e+00 -1.075291e-08 -1.682632e-08 4.340990e-08
## [16,] 0.000000e+00 0.000000e+00 -1.078175e-08 -1.687144e-08 1.682632e-08
## [,16]
## [1,] 0.000000e+00
## [2,] 0.000000e+00
## [3,] 0.000000e+00
## [4,] 0.000000e+00
## [5,] 0.000000e+00
## [6,] 0.000000e+00
## [7,] 0.000000e+00
## [8,] 0.000000e+00
## [9,] 0.000000e+00
## [10,] 0.000000e+00
## [11,] 0.000000e+00
## [12,] 0.000000e+00
## [13,] -1.078175e-08
## [14,] -1.687144e-08
## [15,] 1.682632e-08
## [16,] 1.687144e-08#Calcular de matriz D-lineal
Di <- diag(1/as.vector(pi)) # Contrastar con formula pg 30, Di se plantea de está forma porque para la matriz de varianzas y covarianzas (var.cov.f) se necesita es la inversa.
Di## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 114.9425 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [2,] 0.0000 769.2308 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [3,] 0.0000 0.0000 555.5556 0.000000 0.00 0.000 0.0000 0.000000
## [4,] 0.0000 0.0000 0.0000 1.011941 0.00 0.000 0.0000 0.000000
## [5,] 0.0000 0.0000 0.0000 0.000000 781.25 0.000 0.0000 0.000000
## [6,] 0.0000 0.0000 0.0000 0.000000 0.00 5263.158 0.0000 0.000000
## [7,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 558.6592 0.000000
## [8,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 1.003271
## [9,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [10,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [11,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [12,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [13,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [14,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [15,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [16,] 0.0000 0.0000 0.0000 0.000000 0.00 0.000 0.0000 0.000000
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## [1,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [2,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [3,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [4,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [5,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [6,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [7,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [8,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [9,] 133.3333 0.0000 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [10,] 0.0000 116.2791 0.0000 0.000000 0.000 0.000 0.000 0.00000
## [11,] 0.0000 0.0000 86.2069 0.000000 0.000 0.000 0.000 0.00000
## [12,] 0.0000 0.0000 0.0000 1.028489 0.000 0.000 0.000 0.00000
## [13,] 0.0000 0.0000 0.0000 0.000000 1149.425 0.000 0.000 0.00000
## [14,] 0.0000 0.0000 0.0000 0.000000 0.000 7692.308 0.000 0.00000
## [15,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 5555.556 0.00000
## [16,] 0.0000 0.0000 0.0000 0.000000 0.000 0.000 0.000 10.11941# Calculo de la matriz de varianzas y covarianzas Sigma f gorro
var.cov.ln <- k%*%Di%*%var.cov.lineal%*%Di%*%t(k)
var.cov.ln## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.4448854 0.3667372 0.00000000 0.00000000 0.000000000 0.000000000
## [2,] 0.3667372 0.3587899 0.00000000 0.00000000 0.000000000 0.000000000
## [3,] 0.0000000 0.0000000 0.30929771 0.03744581 0.000000000 0.000000000
## [4,] 0.0000000 0.0000000 0.03744581 0.01012489 0.000000000 0.000000000
## [5,] 0.0000000 0.0000000 0.00000000 0.00000000 0.011769988 0.007221544
## [6,] 0.0000000 0.0000000 0.00000000 0.00000000 0.007221544 0.019739480
## [7,] 0.0000000 0.0000000 0.00000000 0.00000000 0.000000000 0.000000000
## [8,] 0.0000000 0.0000000 0.00000000 0.00000000 0.000000000 0.000000000
## [,7] [,8]
## [1,] 0.0000000 0.0000000
## [2,] 0.0000000 0.0000000
## [3,] 0.0000000 0.0000000
## [4,] 0.0000000 0.0000000
## [5,] 0.0000000 0.0000000
## [6,] 0.0000000 0.0000000
## [7,] 0.8219046 0.6492206
## [8,] 0.6492206 1.3379217# Calculo de la matriz de varianzas y covarianzas Sigma exp f gorro
diagonal <- function(xx){
k<-length(xx)
temp<-matrix(0,ncol=k,nrow=k)
for(i in 1:k)temp[i,i]<-xx[i]
return(temp)
}
Dln <- diagonal(exp(kln))
Dln## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0.87 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00 0.0000000
## [2,] 0.00 0.9981818 0.0000000 0.0000000 0.0000000 0.0000000 0.00 0.0000000
## [3,] 0.00 0.0000000 0.8707483 0.0000000 0.0000000 0.0000000 0.00 0.0000000
## [4,] 0.00 0.0000000 0.0000000 0.9982074 0.0000000 0.0000000 0.00 0.0000000
## [5,] 0.00 0.0000000 0.0000000 0.0000000 0.4658385 0.0000000 0.00 0.0000000
## [6,] 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.9882102 0.00 0.0000000
## [7,] 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.87 0.0000000
## [8,] 0.00 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00 0.9981818var.cov.f.exp <- Q%*%Dln%*%var.cov.ln%*%Dln%*%t(Q)
var.cov.f.exp## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.3367337 0.3184813 0.00000000 0.00000000 0.000000000 0.000000000
## [2,] 0.3184813 0.3574864 0.00000000 0.00000000 0.000000000 0.000000000
## [3,] 0.0000000 0.0000000 0.23451033 0.03254742 0.000000000 0.000000000
## [4,] 0.0000000 0.0000000 0.03254742 0.01008862 0.000000000 0.000000000
## [5,] 0.0000000 0.0000000 0.00000000 0.00000000 0.002554152 0.003324411
## [6,] 0.0000000 0.0000000 0.00000000 0.00000000 0.003324411 0.019276774
## [7,] 0.0000000 0.0000000 0.00000000 0.00000000 0.000000000 0.000000000
## [8,] 0.0000000 0.0000000 0.00000000 0.00000000 0.000000000 0.000000000
## [,7] [,8]
## [1,] 0.0000000 0.000000
## [2,] 0.0000000 0.000000
## [3,] 0.0000000 0.000000
## [4,] 0.0000000 0.000000
## [5,] 0.0000000 0.000000
## [6,] 0.0000000 0.000000
## [7,] 0.6220996 0.563795
## [8,] 0.5637950 1.333061# Matriz de correlaciones
cov2cor(var.cov.f.exp )## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1.0000000 0.9179327 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [2,] 0.9179327 1.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [3,] 0.0000000 0.0000000 1.0000000 0.6691443 0.0000000 0.0000000 0.0000000
## [4,] 0.0000000 0.0000000 0.6691443 1.0000000 0.0000000 0.0000000 0.0000000
## [5,] 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 0.4737772 0.0000000
## [6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.4737772 1.0000000 0.0000000
## [7,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000
## [8,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.6191078
## [,8]
## [1,] 0.0000000
## [2,] 0.0000000
## [3,] 0.0000000
## [4,] 0.0000000
## [5,] 0.0000000
## [6,] 0.0000000
## [7,] 0.6191078
## [8,] 1.0000000library(MASS)
temp <- solve(t(X)%*%ginv(var.cov.f.exp)%*%X)
beta <- temp%*%t(X)%*%ginv(var.cov.f.exp)%*%f
beta## [,1]
## [1,] 0.7741164
## [2,] -0.3078190
## [3,] 0.2193015
## [4,] 0.1499975
## [5,] 0.3680973
## [6,] -0.1298284\(\hat { f^* } =X\hat{\beta}+\epsilon\)
f.gor <- X%*%beta
f.gor ## [,1]
## [1,] 0.7741164
## [2,] 0.9241139
## [3,] 0.9934180
## [4,] 1.0135871
## [5,] 0.4662974
## [6,] 0.9843922
## [7,] 0.6855990
## [8,] 1.0738654