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)
}nro.vars<-2
Nro.Estratos<-12 # Combinacion covariables con P+ y P-
niveles<-c(2,3) # #niveles primer covariable, # niveles segunda covariableGeneración de datos
set.seed(21)
datos <- replicate(6,sample(15:35,4, replace = TRUE))
N <- matrix(datos, 12, 2, byrow = T)###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))
}
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)logistica.bivariable<-function(N,Nro.Estratos,A,X,nro.vars,niveles, imprimir.inter=T){ #Le quite alpha
A<-kronecker(identidad(Nro.Estratos/2),A)
print('Matriz A')
print(A)
#Convertir la matriz N en un vector
N<-matrix(t(N),ncol=1,byrow=F)
print('Tabla de contingencia como vector')
print(N)
# Convertir los ceros que puedan haber en la matriz N a 0.5
# usando la funcion auxiliar ceros.a.5
N<-ceros.a.5(N)
print('ceros a 0.5')
print('N')
#calcular el vector pi-gorro
donde.voy<-0
for(i in 1:Nro.Estratos){
# i=1
identifica<-c('Estrato ',i)
print(identifica)
temp<-N[(donde.voy+1):(donde.voy+2)]
print('Frecuencias de la subpoblación')
print(temp)
donde.voy<-donde.voy+2
n.temp<-sum(temp)
probab<-matrix(temp/n.temp,ncol=1)
print('p estimado del estrato')
print(t(probab))
temp1<- t(probab)
temp2<-rev(probab)
varcov.p<-(matrix(c(temp1, temp2), 2, byrow = T))/n.temp
print('Matriz de varianzas y covarianzas estimada del estrato')
print(varcov.p)
#Matriz de varianza y covarianzas pi-gorro
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)
print('vector que acumula las probabilidades de todos los estratos anteriores')
print(prob.grande)
}
}
#Calcular el vector f = A*pi-gorro
f<-A%*%prob.grande
print('vector f= A*pi-gorro')
print(f)
# Calculo de la matriz de varianzas y covarianzas Sigma f gorro (A*sigma pi-gorro*A')
var.cov.f<-A%*%varcov.grande%*%t(A)
print('var.cov.f gorro')
print(var.cov.f)
# Calcular los betas utilizando inversa generalizada y
library(MASS)
temp<-solve(t(X)%*%ginv(var.cov.f)%*%X)
beta<-temp%*%t(X)%*%ginv(var.cov.f)%*%f
print('Betas estimados')
print(beta)
# matriz de varianzas y covarianzas del modelo
#o matriz de varianzas y covarianzas de beta
Sigma.beta<-temp
print('Matriz de varianzas y covarianzas de Beta')
print(Sigma.beta)
# Calcular estimacion de f-gorro, osea f*
f.gor<-X%*%beta
print('f.gorro*')
print(f.gor)
# error, Z score, p-valor
Z<-nrow(beta)/2
LM<-function(f,X,Sf){
require(MASS)
temp<-t(X)%*%ginv(Sf)
varcov.b<-solve(temp%*%X)
beta<-varcov.b%*%temp%*%f
valores.z<-beta/sqrt(diag(varcov.b))
estadisticos<-cbind(beta,sqrt(diag(varcov.b)),valores.z,2*pnorm(abs(valores.z),lower.tail=F))
colnames(estadisticos)<-c('beta','error','z-valor','p-valor')
estadisticos_sensibilidad <- estadisticos[1:Z,]
print('********** Estimacion del modelo para la sensibilidad *****************')
print(estadisticos_sensibilidad)
print('*********************************************')
estadisticos_especificidad <- estadisticos[(Z+1):nrow(beta),]
print('********** Estimacion del modelo para la especificidad*****************')
print(estadisticos_especificidad)
print('*********************************************')
}
LM(f,X,var.cov.f)
################## Inferencia sobre el modelo
# Pruebas globales
Z<-nrow(beta)/2
#Prueba Global para f1
betaf1<-beta[c(2:Z),]
Sigma.betaf1<-Sigma.beta[c(2:Z),c(2:Z)]
mi.chitempf1g<-t(betaf1)%*%solve(Sigma.betaf1)%*%(betaf1)
#print('Chi cuadrado global para f1')
#print(mi.chitempf1g)
y <- qr(Sigma.betaf1)
gltempgf1g<-y$rank
valor.ptempf1g<-pchisq(mi.chitempf1g,gltempgf1g,lower.tail=F)
#print('Valor p prueba global f1')
#print(valor.ptempf1g)
#Prueba global para f2
P<-Z+1
betaf2<-beta[-c(1,P),]
Sigma.betaf2<-Sigma.beta[-c(1,P),-c(1,P)]
mi.chitempf2g<-t(betaf2)%*%solve(Sigma.betaf2)%*%(betaf2)
#print('Chi cuadrado global para f2')
#print(mi.chitempf2g)
y2 <- qr(Sigma.betaf2)
gltempgf2g<-y2$rank
valor.ptempf2g<-pchisq(mi.chitempf2g,gltempgf2g,lower.tail=F)
#print('Valor p prueba global f2')
#print(valor.ptempf2g)
Tabla_pGlobal <- cbind(c(mi.chitempf1g, mi.chitempf2g), c(valor.ptempf1g, valor.ptempf2g))
colnames(Tabla_pGlobal)<-c('chi-cuadrado','p-valor')
rownames(Tabla_pGlobal)<-c('Sensibilidad','especificidad')
print('********** Pruebas globales *****************')
print(Tabla_pGlobal)
#### Pruebas de coeficientes individuales
#Pruebas chi cuadrado para los parametros individuales para sensibilidad
print('********** Pruebas de coeficientes individuales *****************')
chi.cua<-beta^2/diag(Sigma.beta)
print('Pruebas individuales para beta')
print(chi.cua)
Z <- nrow(beta)/2
print('Z')
print(Z)
chi.cuaf1 <- chi.cua[c(1:Z)]
print('*********************************************')
print('Pruebas individuales para parametros de la sensibilidad')
print(chi.cuaf1)
valor.pf1 <- pchisq(chi.cuaf1,1,lower=F)
#valor.p<-pchisq(chi.cua,1,lower=F) # se redicieron estas dos lineas en la anterior
#valor.pf1<-valor.p[c(1:Z)]
print('Valor-p sensibilidad')
print(valor.pf1)
#Pruebas chi cuadrado para los parametros individuales para especificidad
print('*********************************************')
print('Pruebas individuales para parametros de la especificidad')
Cf2<-identidad(ncol(X)/2)
Cf2<-cbind(Cf2,Cf2)
# print(Sigma.beta) No la imprimi
for(i in 2:nrow(Cf2)){
tempf2 <- matrix(Cf2[i,],nrow=1)
S.temp <- tempf2%*%Sigma.beta%*%t(tempf2)
mi.chitempf2 <- (tempf2%*%beta)%*%solve(S.temp)%*%(tempf2%*%beta)
#gltempf2 <- nrow(tempf2)
#valor.ptempf2 <- pchisq(mi.chitempf2,gltempf2,lower.tail=F)
valor.ptempf2 <- pchisq(mi.chitempf2,1,lower.tail=F) # C es vector por eso el 1 =g= rango de C
# print(paste('Coeficiente', (i-1)))
#print('Coeficiente...')
#print(i-1)
print(paste('Matriz de Pruebas individuales sobre las variables de la especificidad, Coeficiente', (i-1)))
print(tempf2)
print('Chi.obs.....Grados de libertad....Valor-p')
print(c(mi.chitempf2,1,valor.ptempf2))
}
### Pruebas individuales sobre las variables
print("", quote = FALSE)
print('********** Pruebas individuales sobre las variables *****************')
print("", quote = FALSE)
#Pruebas individuales para las variables sensibilidad
for(i in 1:nro.vars){
filas.cv <- niveles[i]-1
Cv <- matrix(0,ncol=ncol(X),nrow=filas.cv)
finalizav <- 1 + sum(niveles[1:i]-1)
iniciav <- finalizav - (niveles[i]-1) + 1
for(j in 1:nrow(Cv)) {
Cv[j,iniciav-1+j] <- 1
}
# mi.chiv<-t(Cv%*%beta)%*%solve(Cv%*%Sigma.beta%*%t(Cv))%*%Cv%*%beta
mi.chiv <- t(Cv%*%beta)%*%solve(Cv%*%Sigma.beta%*%t(Cv))%*%Cv%*%beta
gl <- nrow(Cv)
valor.p <- pchisq(mi.chiv,gl,lower.tail=F)
#print('Variable...')
#print(i)
#print('Matriz de Pruebas individuales sobre las variables de la sensibilidad')
print(paste('Matriz de Pruebas individuales de la sensibilidad con respecto a la variable', i))
print(Cv)
print('Chi.obs.....Grados de libertad....Valor-p')
print(c(mi.chiv,gl,valor.p))
} # for i
#Pruebas individuales para las variables especificidad
print("", quote = FALSE)
for(i in 1:nro.vars){
filas.cv <- niveles[i] - 1
Cv1 <- matrix(0,ncol=ncol(X)/2,nrow=filas.cv)
finalizav1 <- 1 + sum(niveles[1:i]-1)
iniciav1 <- finalizav1 - (niveles[i]-1) + 1
for(j in 1:nrow(Cv1)) {Cv1[j,iniciav1-1+j] <- 1}
Cv2 <- matrix(0,ncol=ncol(X)/2,nrow=filas.cv)
finalizav2 <- 1 + sum(niveles[1:i]-1)
iniciav2 <- finalizav2 - (niveles[i]-1) + 1
for(j in 1:nrow(Cv2)) {Cv2[j,iniciav1-1+j] <- 1}
Cvf2 <- cbind(Cv1,Cv2)
mi.chicvf2 <- t(Cvf2%*%beta)%*%solve(Cvf2%*%Sigma.beta%*%t(Cvf2))%*%Cvf2%*%beta
glcvf2 <- nrow(Cvf2)
valor.pcvf2 <- pchisq(mi.chicvf2,glcvf2,lower.tail=F)
#print('Variable...')
#print(i)
print(paste('Matriz de Pruebas individuales de la especificidad sobre la variable ', i))
print(Cvf2)
print('Chi.obs.....Grados de libertad....Valor-p')
print(c(mi.chicvf2,glcvf2,valor.pcvf2))
} # for i
## Pruebas de Contraste sobre las Variables
print("", quote = FALSE)
print('********** Pruebas de Contraste sobre las Variables *****************')
# Contrastes para la sensibilidad
if(any(niveles>2)){
for(i in 1:nro.vars){
if(niveles[i]>2){
filas.con<-niveles[i]-2
# Si la matriz de dise?o no tiene incluido el intercepto use...
# Con<-matrix(0,ncol=ncol(X),nrow=filas.con)
# Si la matriz de dise?o tiene incluido el intercepto use...
Con<-matrix(0,ncol=ncol(X),nrow=filas.con)
finaliza <- 1+sum(niveles[1:i]-1)
inicia <- finaliza-(niveles[i]-1)+1
Con[1:nrow(Con),inicia] <- 1
for(j in 1:nrow(Con)) {Con[j,inicia+j] <- -1}
mi.chi <- t(Con%*%beta)%*%solve(Con%*%Sigma.beta%*%t(Con))%*%Con%*%beta
gl <- nrow(Con)
valor.p <- pchisq(mi.chi,gl,lower.tail=F)
#print(paste('Variable...',i))
#print(i)
#print('Matriz de contrastes')
print("", quote = FALSE)
print("Sensibilidad", quote = FALSE)
print(paste('Matriz de contrastes variabls', i))
print(Con)
print('Chi.obs.....Grados de libertad....Valor-p')
print(c(mi.chi,gl,valor.p))
} # Fin if
} # for i
} #fin if any
# Contrastes para la especificidad
if(any(niveles>2)){
for(i in 1:nro.vars){
if(niveles[i]>2){
filas.con <- niveles[i]-2
Con1 <- matrix(0,ncol=ncol(X)/2,nrow=filas.con)
finaliza <- 1+sum(niveles[1:i]-1)
inicia <- finaliza-(niveles[i]-1)+1
Con1[1:nrow(Con1),inicia] <- 1
for(j in 1:nrow(Con1))Con1[j,inicia+j] <- -1
Con2 <- matrix(0,ncol=ncol(X)/2,nrow=filas.con)
finaliza2 <- 1+sum(niveles[1:i]-1)
inicia2 <- finaliza2-(niveles[i]-1)+1
Con2[1:nrow(Con2),inicia2] <- 1
for(j in 1:nrow(Con2)) {Con2[j,inicia+j]<- -1}
Conf2 <- cbind(Con1,Con2)
mi.chi2 <- t(Conf2%*%beta)%*%solve(Conf2%*%Sigma.beta%*%t(Conf2))%*%Conf2%*%beta
gl2 <- nrow(Conf2)
valor.p2 <- pchisq(mi.chi2,gl2,lower.tail=F)
#print(paste('Variable...',i))
#print(i)
print("", quote = FALSE)
print("Especificidad", quote = FALSE)
print(paste('Matriz de contrastes variable', i))
print(Conf2)
print('Chi.obs.....Grados de libertad....Valor-p')
print(c(mi.chi2,gl2,valor.p2))
} # Fin if
} # for i
} #fin if any
}# Fin funcion principallogistica.bivariable(N=N,Nro.Estratos=12,A=A,X=X,nro.vars=2,niveles=c(2,3))