Matrices a usar internamente

La función crea.bloque permite crear una matriz en bloque a partir de dos matrices:

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)
}

La funcion ceros.a.5 corrige los ceros en 0.5

ceros.a.5<-function(Tabla){
  Nro.col<-ncol(Tabla)
  if(any(Tabla==0))Tabla<-matrix(ifelse
                                 (as.vector(Tabla)==0,0.5,Tabla),ncol=Nro.col)
  return(Tabla)
}

La funcion identidad permite construir una matriz diagonal con k elementos iguales a 1 en su diagonal principal

identidad<-function(k)diag(rep(1,k))

La funcion bloque permite construir una matriz donde su primera fila esta compuesta de ceros y debajo se adiciona una matriz diagonal con k elementos iguales a 1 en su diagonal principal.

bloque<-function(k)rbind(0,identidad(k))

La funcion repita toma dos matrices apartir de las cuales se genera una nueva matriz de mayor dimension. La primera matriz aparece en las primeras columnas y la segunda se adiciona ocupando las ultimas columnas de la nueva matriz generada. Cada fila aparece repetida k veces una debajo de la otra.

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)
}

GSK

Matriz de Diseño

matriz.diseño<-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))
}

##################### Hablar de estos datos en general

X <- matriz.diseño(2,c(2,3))
set.seed(21)
datos <- replicate(6,sample(15:35,4, replace = TRUE)) 
N <- matrix(datos, 12, 2, byrow = T)

nro.vars<-2

A <- matrix(c(1,0,0,0,
              0,0,0,1),ncol=4,byrow=T)

Nro.Estratos<-12
niveles<-c(2,3)

Funciones propuestas

1. De manera directa

####  1. construccion del Vector de probabilidades por subpoblacion $\pi$
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)

# Vector con las probabilidades de las tablas de todas las subpoblaciones 
pi <- c(pi_pob1, pi_pob2, pi_pob3, pi_pob4)




#### 2. Constrir vector $A\pi$

# Matriz A bloque interno para un estrato

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


# Matriz A completa para todos los estratos
A <- kronecker(identidad(Nro.Estratos),A) 
A


# Vector A*pi
A_Pi <- A%*%pi
A_Pi


#### 3. Vector ln(A*pi)

LogNatu <- log(A_Pi)
LogNatu


#### 4. Vector k*ln(A*pi)

# Construccion matriz K
k <- matrix(c(1, -1, 0, 0,
              0, 0, 1, -1),2,4, byrow = T)

k <- kronecker(identidad(Nro.Estratos),k) 
k

# Vector k*ln(A*pi)
kln <- k%*%LogNatu
kln


#### 5. Vector exp(k*ln(A*pi))

ekln <- exp(kln)
ekln

#### 6. Vector Q*exp(k*ln(A*pi))

# Constricción matriz Q
Q <- diag(1,2*Nro.Estratos)
Q

# Calcular f_exp-estimado = Q*exp(k*ln(A*pi))
f <- Q%*%ekln
f


################################ Estimación de los parámetros
# Para estimar $\beta$, comenzamos a hallar o definir:

#### 1. Construir matriz de diseño para el Modelo General X = OoX

X <- matriz.diseno(nro.vars, niveles)


#### 2. Matriz de varianzas y covarianzas ${ \hat { \Sigma } }_{ \hat { f } }$


N <- pob

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
  }
}


# Calcular (A*sigma pi-gorro*A")
var.cov.lineal <- A%*%varcov.grande%*%t(A)
var.cov.lineal

#Calcular de matriz D-lineal
Di <- diag(1/as.vector(pi))           

# 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

# 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

var.cov.f.exp <- Q%*%Dln%*%var.cov.ln%*%Dln%*%t(Q)
var.cov.f.exp

# Matriz de correlaciones
cov2cor(var.cov.f.exp )



#### 3. Estimación de los $\beta$ 

library(MASS)

temp <- solve(t(X)%*%ginv(var.cov.f.exp)%*%X)
beta <- temp%*%t(X)%*%ginv(var.cov.f.exp)%*%f
beta


#### 4. Estimación de la sensibilidad y especificidad por subpoblaciones

# $f^{*}$ estimado, $\hat { f^* } =X\hat{\beta}+\epsilon$

f.gor <- X%*%beta
f.gor