Aquí se definen parámetros generales y el espacio de trabajo
options(max.print = 100000, width=300) # dimensiones salida
getwd() # donde estamos ahora
## [1] "/home/sergio/Documents/Lavoro in corso/Análisis de datos y estadística inferencial - IMCA/Lezioni 2024/Clase_5/R"
#setwd("~/Documents/Lavoro in corso/Análisis de datos y estadística inferencial - IMCA/Lezioni 2024/Clase_5 /R") # donde vamos
#Cruzado de variables cualitativas
##Variables dicotómicas
Si se tienen dos variables dicotómicas, se construye una tabla como la siguiente, que cruza el Infarto con el uso de aspirina o de placebo (estudio longitudinal). Si ya existe, se puede volver a los datos brutos como sigue:
# Infarto No_infarto
# Aspirina 104 10933
# Placebo 189 10845
################################
a <- 104
b <- 10933
c <- 189
d <- 10845
n <- a+b+c+d
data <- matrix(rep(0,n*2),ncol=2,nrow=n)
for (i in 1:a) {data[i,] <- c("Aspirina","Infarto")}
for (i in (a+1):(a+b)) {data[i,] <- c("Aspirina","No_infarto")}
for (i in (a+b+1):(a+b+c)) {data[i,] <- c("Placebo","Infarto")}
for (i in (a+b+c+1):n) {data[i,] <- c("Placebo","No_infarto")}
Así ahora suponemos de tener los datos brutos.
t <- table(data[,1],data[,2]) # tabla
noquote("Tabla")
## [1] Tabla
t
##
## Infarto No_infarto
## Aspirina 104 10933
## Placebo 189 10845
tm <- addmargins(t) # más margenes
noquote("Tabla con márgenes")
## [1] Tabla con márgenes
tm
##
## Infarto No_infarto Sum
## Aspirina 104 10933 11037
## Placebo 189 10845 11034
## Sum 293 21778 22071
f <- t/n # tabla de frecuencias
fm <- addmargins(f) # más márgenes
noquote("Tabla de frecuencias")
## [1] Tabla de frecuencias
round(fm,3)
##
## Infarto No_infarto Sum
## Aspirina 0.005 0.495 0.500
## Placebo 0.009 0.491 0.500
## Sum 0.013 0.987 1.000
riesgo <- fm[3,1]/fm[3,3] # riesgo
noquote(paste0("Riesgo = ",round(riesgo,5)))
## [1] Riesgo = 0.01328
Odd <- fm[1,1]/fm[1,2] # momio
noquote(paste0("Momio = ",round(Odd,5)))
## [1] Momio = 0.00951
El riesgo de infarto 0.013/0.987 es muy bajo. Los riesgos condicionados por Aspirina/ Placebo son:
rI_A <- fm[1,1]/fm[1,3] # riesgo con aspirina
rI_P <- fm[2,1]/fm[2,3] # riesgo sin aspirina
riesgos <- c(rI_A,rI_P)
names(riesgos) <- colnames(t)
noquote("Riesgos condicionados")
## [1] Riesgos condicionados
round(riesgos,5) # riesgos condicionados
## Infarto No_infarto
## 0.00942 0.01713
RR <- rI_A/rI_P # razón de riesgo
noquote(paste0("Razón de riesgos = ",round(RR,5)))
## [1] Razón de riesgos = 0.55011
OddI_A <- fm[1,1]/fm[1,2] # momio con aspirina
OddI_P <- fm[2,1]/fm[2,2] # momio sin aspirina
momios <- c(OddI_A,OddI_P)
names(momios) <- colnames(t)
noquote("Momios condicionados")
## [1] Momios condicionados
round(momios,5)
## Infarto No_infarto
## 0.00951 0.01743
OR <- OddI_A/OddI_P # razón de momios
noquote(paste0("Razón de momios = ",round(OR,5)))
## [1] Razón de momios = 0.54584
Se resulta que el riesgo de infarto con aspirina se reduce de cuasi la mitad. La razón de momios se exprime como alrededor de un infarto con aspirina contra dos con placebo. En ambos casos se resulta que la aspirina hace bien. En este caso, siendo un estudio longitudinal, se puede hacer referencia a los riesgos.
Se puede hacer una función que calcula riesgos y momios. Notese que hay que darse una regla de como está formada la tabla, como sigue:\
factor de riesgo
presente ausente
tratamiento
no tratamiento
risk_odds <- function(riesgo,tratamiento) {
res <- list() # list for results
# tables
t <- table(tratamiento,riesgo) # contingency table
tm <- addmargins(t) # with margins
res$Tabla <- tm # table with margins
f <- t/n # frequencies
fm <- addmargins(f) # with margins
res$Frecuencias <- round(fm,3) # table with margins
# risks
risk <- fm[3,1]/fm[3,3] # factor risk
rI_T <- fm[1,1]/fm[1,3] # risk for treatment
rI_N <- fm[2,1]/fm[2,3] # risk fro no treatment
RR <- rI_T/rI_N # risks ratio
tabrisk <- as.table(matrix(c(risk,NA,rI_T,rI_N,RR,NA),nrow=3,ncol=2,byrow=TRUE))
rownames(tabrisk) <- c("riesgo","condicionados","razón de riesgos")
colnames(tabrisk) <- rownames(t)
res$Riesgos <- print(round(tabrisk,5), na.print = "" , quote = FALSE)
# odds
Odd <- fm[1,1]/fm[1,2] # odd
OddI_T <- fm[1,1]/fm[1,2] # odd for treatment
OddI_N <- fm[2,1]/fm[2,2] # odd for no treatment
OR <- OddI_T/OddI_N # odds ratio
tabodds <- as.table(matrix(c(Odd,NA,OddI_T,OddI_N,OR,NA),nrow=3,ncol=2,byrow=TRUE))
rownames(tabodds) <- c("momio","condicionados","razón de momios")
colnames(tabodds) <- rownames(t)
res$Momios <- print(round(tabodds,5), na.print = "" , quote = FALSE)
res
}
Así, corriendo los datos de Infarto, se resulta
risk_odds(data[,2],data[,1])
## Aspirina Placebo
## riesgo 0.01328
## condicionados 0.00942 0.01713
## razón de riesgos 0.55011
## Aspirina Placebo
## momio 0.00951
## condicionados 0.00951 0.01743
## razón de momios 0.54584
## $Tabla
## riesgo
## tratamiento Infarto No_infarto Sum
## Aspirina 104 10933 11037
## Placebo 189 10845 11034
## Sum 293 21778 22071
##
## $Frecuencias
## riesgo
## tratamiento Infarto No_infarto Sum
## Aspirina 0.005 0.495 0.500
## Placebo 0.009 0.491 0.500
## Sum 0.013 0.987 1.000
##
## $Riesgos
## Aspirina Placebo
## riesgo 0.01328
## condicionados 0.00942 0.01713
## razón de riesgos 0.55011
##
## $Momios
## Aspirina Placebo
## momio 0.00951
## condicionados 0.00951 0.01743
## razón de momios 0.54584
Consideramos la tabla de Snee, que cruza color de ojos y color del pelo de 592 alumnos. Transformamos los datos en factores ordenados, así que se guarda un orden, evitando que se pongan por orden alfabético.
snee=read.csv("Snee.csv") # se trabaja con datos de Snee
snee$Color_ojos <- factor(snee$Color_ojos, # factor y etiquetas de la primera variable
levels = c(1,2,3,4),
labels = c("Oscuros","Pardos",
"Verdes","Azules"),ordered=TRUE)
snee$Color_pelo <- factor(snee$Color_pelo, # factor y etiquetas de la segunda variable
levels = c(1,2,3,4),
labels = c("Negro","Castano",
"Rojo","Rubio"),ordered=TRUE)
Se construye una tabla de contingencia con márgenes, y de frecuencias
# datos
tsnee <- table(snee[,1],snee[,2]); tsnee
##
## Negro Castano Rojo Rubio
## Oscuros 68 119 26 7
## Pardos 15 54 14 10
## Verdes 5 29 14 16
## Azules 20 84 17 94
nr <- dim(tsnee)[1] # número de filas
nc <- dim(tsnee)[2] # número de columnas
totsnee <- sum(tsnee); totsnee # total de unidades
## [1] 592
# totales marginales
tr <- margin.table(tsnee,1);tr # totales marginales de fila
##
## Oscuros Pardos Verdes Azules
## 220 93 64 215
tc <- margin.table(tsnee,2);tc # y de columna
##
## Negro Castano Rojo Rubio
## 108 286 71 127
tcsnee <- addmargins(tsnee); tcsnee # se van ajuntanto los totales
##
## Negro Castano Rojo Rubio Sum
## Oscuros 68 119 26 7 220
## Pardos 15 54 14 10 93
## Verdes 5 29 14 16 64
## Azules 20 84 17 94 215
## Sum 108 286 71 127 592
colnames(tcsnee)[nc+1] <- "Tot_ojos" # con su nombres
rownames(tcsnee)[nr+1] <- "Tot_pelo"
# tabla completa
print("Tabla de contingencia") # y se imprime
## [1] "Tabla de contingencia"
tcsnee
##
## Negro Castano Rojo Rubio Tot_ojos
## Oscuros 68 119 26 7 220
## Pardos 15 54 14 10 93
## Verdes 5 29 14 16 64
## Azules 20 84 17 94 215
## Tot_pelo 108 286 71 127 592
str(tcsnee)
## 'table' num [1:5, 1:5] 68 15 5 20 108 119 54 29 84 286 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:5] "Oscuros" "Pardos" "Verdes" "Azules" ...
## ..$ : chr [1:5] "Negro" "Castano" "Rojo" "Rubio" ...
# tabla de frecuencias
ptsnee <- tcsnee/totsnee # tabla de frecuencia con márgenes
print("Tabla de frecuencias")
## [1] "Tabla de frecuencias"
round(ptsnee, 3)
##
## Negro Castano Rojo Rubio Tot_ojos
## Oscuros 0.115 0.201 0.044 0.012 0.372
## Pardos 0.025 0.091 0.024 0.017 0.157
## Verdes 0.008 0.049 0.024 0.027 0.108
## Azules 0.034 0.142 0.029 0.159 0.363
## Tot_pelo 0.182 0.483 0.120 0.215 1.000
# tabla sin márgenes
psnee <- ptsnee[c(1:nr),c(1:nc)] # sin margenes
pr <- ptsnee[1:nr,nc+1]; pr # margenes
## Oscuros Pardos Verdes Azules
## 0.3716216 0.1570946 0.1081081 0.3631757
pc <- ptsnee[nr+1,1:nc]; pc
## Negro Castano Rojo Rubio
## 0.1824324 0.4831081 0.1199324 0.2145270
round(psnee,3)
##
## Negro Castano Rojo Rubio
## Oscuros 0.115 0.201 0.044 0.012
## Pardos 0.025 0.091 0.024 0.017
## Verdes 0.008 0.049 0.024 0.027
## Azules 0.034 0.142 0.029 0.159
Aquí se calculan las distribuciones de cada variable condicionadas para otra.
# perfiles de filas y su gráfico
prsnee <- prop.table(tsnee,1) # perfíl de fila
print("Perfiles de filas")
## [1] "Perfiles de filas"
round(prsnee,3)
##
## Negro Castano Rojo Rubio
## Oscuros 0.309 0.541 0.118 0.032
## Pardos 0.161 0.581 0.151 0.108
## Verdes 0.078 0.453 0.219 0.250
## Azules 0.093 0.391 0.079 0.437
apply(prsnee,1,sum) # control
## Oscuros Pardos Verdes Azules
## 1 1 1 1
prfr <- rbind(prsnee,pc) # impresión
round(cbind(prfr,apply(prfr,1,sum)),3)
## Negro Castano Rojo Rubio
## Oscuros 0.309 0.541 0.118 0.032 1
## Pardos 0.161 0.581 0.151 0.108 1
## Verdes 0.078 0.453 0.219 0.250 1
## Azules 0.093 0.391 0.079 0.437 1
## pc 0.182 0.483 0.120 0.215 1
# grafico
plot(1:nc,ptsnee[nr+1,1:nc],ylim=c(0,max(prsnee+.05)),type="n", xaxt="n", # y su gráfico
xlab = "Color of hair", ylab = "Frequencies",
main = "Snee data: profiles of hair", sub = "conditioned by eyes")
axis(1,at=c(1:nc),labels=colnames(psnee)[1:nc]) # eje horizontal
lines(1:nc,ptsnee[5,1:nc],col="red") # línias de perfiles
lines(1:nc,prsnee[1,1:nc],col="black")
lines(1:nc,prsnee[2,1:nc],col="brown")
lines(1:nc,prsnee[3,1:nc],col="green")
lines(1:nc,prsnee[4,1:nc],col="blue")
legend("topright",legend=c(rownames(psnee),"Marginal"),lty=c(1,1),cex=0.6,y.intersp=.7,
col=c("black","brown","green","blue","red")) # leyenda
# perfiles de columnas y su gráfico
pcsnee <- prop.table(tsnee,2) # perfil de columnas
print("Perfiles de columnas")
## [1] "Perfiles de columnas"
round(pcsnee,3)
##
## Negro Castano Rojo Rubio
## Oscuros 0.630 0.416 0.366 0.055
## Pardos 0.139 0.189 0.197 0.079
## Verdes 0.046 0.101 0.197 0.126
## Azules 0.185 0.294 0.239 0.740
apply(pcsnee,2,sum) # control
## Negro Castano Rojo Rubio
## 1 1 1 1
# grafico
plot(1:nr,ptsnee[nc+1,1:nr],ylim=c(0,max(pcsnee+0.05)),type="n", xaxt="n", # y su gráfico
xlab = "Color of eyes", ylab = "Frequencies",
main = "Snee data: profiles of eyes", sub = "conditioned by hair")
axis(1,at=c(1:nr),labels=rownames(psnee)[1:nr]) # eje horizontal
lines(1:nr,ptsnee[1:nr,nc+1],col="blue") # linias de perfiles
lines(1:nr,pcsnee[1:nr,1],col="black")
lines(1:nr,pcsnee[1:nr,2],col="brown")
lines(1:nr,pcsnee[1:nr,3],col="red")
lines(1:nr,pcsnee[1:nr,4],col="orange")
legend(3,0.8,legend=c(colnames(psnee),"Marginal"),lty=c(1,1),cex=0.6,y.intersp=.7,
col=c("black","brown","red","orange","blue")) # leyenda
# desvíos a los perfiles marginales
print("Desvíos a los perfiles de fila")
## [1] "Desvíos a los perfiles de fila"
round(rbind(sweep(prsnee,2,pc,"-"),pr),3) # desvíos más marginal
## Negro Castano Rojo Rubio
## Oscuros 0.127 0.058 -0.002 -0.183
## Pardos -0.021 0.098 0.031 -0.107
## Verdes -0.104 -0.030 0.099 0.035
## Azules -0.089 -0.092 -0.041 0.223
## pr 0.372 0.157 0.108 0.363
print("Desvíos a los perfiles de columna")
## [1] "Desvíos a los perfiles de columna"
round(cbind(sweep(pcsnee,1,pr,"-"),pc),3) # desvíos más marginal
## Negro Castano Rojo Rubio pc
## Oscuros 0.258 0.044 -0.005 -0.317 0.182
## Pardos -0.018 0.032 0.040 -0.078 0.483
## Verdes -0.062 -0.007 0.089 0.018 0.120
## Azules -0.178 -0.069 -0.124 0.377 0.215
# tabla de frecuencias esperadas
exsnee <- pr%*%t(pc); exsnee # tabla valores esperados
## Negro Castano Rojo Rubio
## [1,] 0.06779584 0.17953342 0.04456949 0.07972288
## [2,] 0.02865915 0.07589367 0.01884074 0.03370104
## [3,] 0.01972243 0.05222790 0.01296567 0.02319211
## [4,] 0.06625502 0.17545311 0.04355654 0.07791100
rownames(exsnee) == rownames(tsnee)
## logical(0)
# los perfiles son iguales
prexsnee <- prop.table(exsnee,1);prexsnee # perfiles de fila
## Negro Castano Rojo Rubio
## [1,] 0.1824324 0.4831081 0.1199324 0.214527
## [2,] 0.1824324 0.4831081 0.1199324 0.214527
## [3,] 0.1824324 0.4831081 0.1199324 0.214527
## [4,] 0.1824324 0.4831081 0.1199324 0.214527
pcexsnee <- prop.table(exsnee,2);pcexsnee # perfiles de columna
## Negro Castano Rojo Rubio
## [1,] 0.3716216 0.3716216 0.3716216 0.3716216
## [2,] 0.1570946 0.1570946 0.1570946 0.1570946
## [3,] 0.1081081 0.1081081 0.1081081 0.1081081
## [4,] 0.3631757 0.3631757 0.3631757 0.3631757
pexsnee <- addmargins(exsnee) # se ajuntan los totales
colnames(pexsnee)[nc+1] <- "Tot_ojos"
rownames(pexsnee)[nr+1] <- "Tot_pelo"
print("Tabla de frecuencias esperadas")
## [1] "Tabla de frecuencias esperadas"
round(pexsnee,3)
## Negro Castano Rojo Rubio Tot_ojos
## 0.068 0.180 0.045 0.080 0.372
## 0.029 0.076 0.019 0.034 0.157
## 0.020 0.052 0.013 0.023 0.108
## 0.066 0.175 0.044 0.078 0.363
## Tot_pelo 0.182 0.483 0.120 0.215 1.000
print("Tabla de desvíos (frecuencias relativas)")
## [1] "Tabla de desvíos (frecuencias relativas)"
despsnee <- (psnee-exsnee); round(despsnee,3) # desvíos a la independencia
##
## Negro Castano Rojo Rubio
## Oscuros 0.047 0.021 -0.001 -0.068
## Pardos -0.003 0.015 0.005 -0.017
## Verdes -0.011 -0.003 0.011 0.004
## Azules -0.032 -0.034 -0.015 0.081
desnee <- despsnee * totsnee
print("Tabla de desvíos")
## [1] "Tabla de desvíos"
round(desnee,3) # tabla de desvios
##
## Negro Castano Rojo Rubio
## Oscuros 27.865 12.716 -0.385 -40.196
## Pardos -1.966 9.071 2.846 -9.951
## Verdes -6.676 -1.919 6.324 2.270
## Azules -19.223 -19.868 -8.785 47.877
sum(desnee)
## [1] 6.661338e-15
# estadísticas
entr <- -sum(pr*log(pr,2)); entr # entropía filas
## [1] 1.827862
entrr <- entr/log(nr,2); entrr # entropía relativa
## [1] 0.9139308
entc <- -sum(pc*log(pc,2)); entc # entropía columnas
## [1] 1.798227
entcr <- entc/log(nc,2); entcr # entropía relativa
## [1] 0.8991135
entt <- -sum(psnee*log(psnee,2)); entt # entropía conjunta
## [1] 3.447648
enttr <- entt/log(nr*nc,2); enttr # entropía relativa
## [1] 0.861912
infm <- sum(psnee*(log(psnee,2)- log(exsnee,2))); infm # información mútua
## [1] 0.1784404
infm <- sum(psnee*(log(psnee/exsnee,2))); infm # información mútua
## [1] 0.1784404
infm <- entr + entc - entt; infm # información mútua
## [1] 0.1784404
G2 <- 2 * totsnee * infm * log(exp(1),2) # g²
entconr <- entr - infm; entconr # entropías condiciondas
## [1] 1.649421
entconc <- entc - infm; entconc
## [1] 1.619787
# tabla de entropias
entropias <- rbind(c(entr,entrr),c(entc,entcr),c(entt,enttr))
rownames(entropias) <- c("filas","columnas","conjunta")
colnames(entropias) <- c("entropía","relativa")
print("Entropías")
## [1] "Entropías"
round(entropias,5)
## entropía relativa
## filas 1.82786 0.91393
## columnas 1.79823 0.89911
## conjunta 3.44765 0.86191
print("Entropías condicionadas")
## [1] "Entropías condicionadas"
econ <- c(entconr,entconc)
names(econ) <- c("fila|columnas","columna|filas")
round(econ,5)
## fila|columnas columna|filas
## 1.64942 1.61979
noquote(paste0("información mutua = ",infm))
## [1] información mutua = 0.17844039224729
noquote(paste0("G² = ",G2))
## [1] G² = 304.803121683505
# estadísticas sobre chi-cuadrado
phi <- sum((psnee-exsnee)^2/exsnee); phi # phi cuadrado
## [1] 0.2335977
chi <- totsnee * phi; chi # chi cuadrado # chi cuadrado
## [1] 138.2898
df <- (nr-1)*(nc-1) # grados de libertad
pearson <- sqrt(chi/(totsnee+chi)); pearson # C de Pearson
## [1] 0.4351585
tchuprow <- sqrt(chi/(totsnee * sqrt(df))); tchuprow # T de Tchuprow
## [1] 0.2790446
cramer <- sqrt(chi/(totsnee * min((nr-1),(nc-1)))); cramer # phi de Cramer
## [1] 0.2790446
# tabla de estadísticas sobre chi-cuadrado
chist <- rbind(phi,chi,df,pearson,tchuprow,cramer)
rownames(chist) <- c("phi-cuadrado","chi-cuadrado","grados de libertad",
"C de Pearson","T de Tchuprow",
"phi de Cramer")
colnames(chist) <- "estadísticas"
print("Estadísticas sobre chi-cuadrado")
## [1] "Estadísticas sobre chi-cuadrado"
chist
## estadísticas
## phi-cuadrado 0.2335977
## chi-cuadrado 138.2898416
## grados de libertad 9.0000000
## C de Pearson 0.4351585
## T de Tchuprow 0.2790446
## phi de Cramer 0.2790446
Por esto también se puede hacer una función
stat_tab <- function(data) {
res <- list()
# datos
tdata <- table(data[,1],data[,2])
nr <- dim(tdata)[1] # número de filas
nc <- dim(tdata)[2] # número de columnas
n <- sum(tdata) # total de unidades
# totales marginales
tr <- margin.table(tdata,1) # totales marginales de fila
tc <- margin.table(tdata,2) # y de columna
tcdata <- addmargins(tdata) # se van ajuntanto los totales
colnames(tcdata)[nc+1] <- "Margin" # con su nombres
rownames(tcdata)[nr+1] <- "Margin"
# tabla completa
res$Table <- tcdata
# tabla de frecuencias
ptdata <- tcdata/n # tabla de frecuencia con márgenes
res$Frequencies <- round(ptdata, 3)
# tabla sin márgenes
pdata <- ptdata[c(1:nr),c(1:nc)] # sin margenes
pr <- ptdata[1:nr,nc+1] # margenes
pc <- ptdata[nr+1,1:nc]
# perfiles de filas
prdata <- prop.table(tdata,1) # perfíl de fila
prfr <- rbind(prdata,pc) # margenes
prfr <- cbind(prfr,apply(prfr,1,sum))
rownames(prfr)[nr+1] <- "Margin"
colnames(prfr)[nc+1] <- "Margin"
res$row.profiles <- round(prfr,3)
# perfiles de columnas
pcdata <- prop.table(tdata,2) # perfil de columnas
prfc <- cbind(pcdata,pr) # márgenes
prfc <- rbind(prfc,apply(prfc,2,sum))
rownames(prfc)[nr+1] <- "Margin"
colnames(prfc)[nc+1] <- "Margin"
res$col.profiles <- round(prfc,3)
# desvíos a los perfiles marginales
# print("Desvíos a los perfiles de fila")
# round(rbind(sweep(prdata,2,pc,"-"),pr),3) # desvíos más marginal
# print("Desvíos a los perfiles de columna")
# round(cbind(sweep(pcdata,1,pr,"-"),pc),3) # desvíos más marginal
# tabla de frecuencias esperadas
exdata <- pr%*%t(pc) # tabla valores esperados
pexdata <- addmargins(exdata) # se ajuntan los totales
colnames(pexdata)[nc+1]<- "Margin"
rownames(pexdata)[nr+1]<- "Margin"
res$exp.frequencies <- round(pexdata,3)
despdata <- (pdata-exdata)
res$dev.frequencies <- round(despdata,3) # desvíos a la independencia
dedata <- despdata * n # números de desvíos
res$dev.values <- round(dedata,3)
# estadísticas
entr <- -sum(pr*log(pr,2)) # entropía filas
entrr <- entr/log(nr,2) # entropía relativa
entc <- -sum(pc*log(pc,2)) # entropía columnas
entcr <- entc/log(nc,2) # entropía relativa
entt <- -sum(pdata*log(pdata,2)) # entropía conjunta
enttr <- entt/log(nr*nc,2) # entropía relativa
entexp <- -sum(exdata*log(exdata,2)) # entropía esperada
entexpr <- entexp/log(nr*nc,2) # entropía relativa
infm <- entr + entc - entt # información mútua
G2 <- 2 * n * infm * log(2) # g²
entconr <- entr - infm # entropias condicionadas
entconc <- entc - infm
# tabla de entropias
entropias <- rbind(c(entr,entrr),c(entc,entcr),c(entt,enttr),c(entexp,entexpr))
rownames(entropias) <- c("filas","columnas","conjunta","esperada")
colnames(entropias) <- c("entropía","relativa")
res$Entropies <- round(entropias,5)
econ <- c(entconr,entconc)
names(econ) <- c("fila|columnas","columna|filas")
res$Cond.entropies <- round(econ,5)
info <- c(infm,G2)
names(info) <- c("mutual_information","G² statistics")
res$Information <- round(info,5)
# estadísticas sobre chi-cuadrado
phi <- sum((pdata-exdata)^2/exdata) # phi cuadrado
chi <- n * phi # chi cuadrado
df <- (nr-1)*(nc-1) # grados de libertad
pearson <- sqrt(chi/(n+chi)) # pearson
tchuprow <- sqrt(chi/(n * sqrt(df))) # tchuprow
cramer <- sqrt(phi/(min((nr-1),(nc-1)))) # cramer
# tabla de estadísticas sobre chi-cuadrado
chist <- rbind(phi,chi,df,pearson,tchuprow,cramer)
rownames(chist) <- c("phi-cuadrado","chi-cuadrado","grados de libertad",
"C de Pearson","T de Tchuprow",
"phi de Cramer")
colnames(chist) <- "estadísticas"
res$Phi_stats <- round(chist,5)
res
}
Intentamos con los datos de Snee:
stat_tab(snee)
## $Table
##
## Negro Castano Rojo Rubio Margin
## Oscuros 68 119 26 7 220
## Pardos 15 54 14 10 93
## Verdes 5 29 14 16 64
## Azules 20 84 17 94 215
## Margin 108 286 71 127 592
##
## $Frequencies
##
## Negro Castano Rojo Rubio Margin
## Oscuros 0.115 0.201 0.044 0.012 0.372
## Pardos 0.025 0.091 0.024 0.017 0.157
## Verdes 0.008 0.049 0.024 0.027 0.108
## Azules 0.034 0.142 0.029 0.159 0.363
## Margin 0.182 0.483 0.120 0.215 1.000
##
## $row.profiles
## Negro Castano Rojo Rubio Margin
## Oscuros 0.309 0.541 0.118 0.032 1
## Pardos 0.161 0.581 0.151 0.108 1
## Verdes 0.078 0.453 0.219 0.250 1
## Azules 0.093 0.391 0.079 0.437 1
## Margin 0.182 0.483 0.120 0.215 1
##
## $col.profiles
## Negro Castano Rojo Rubio Margin
## Oscuros 0.630 0.416 0.366 0.055 0.372
## Pardos 0.139 0.189 0.197 0.079 0.157
## Verdes 0.046 0.101 0.197 0.126 0.108
## Azules 0.185 0.294 0.239 0.740 0.363
## Margin 1.000 1.000 1.000 1.000 1.000
##
## $exp.frequencies
## Negro Castano Rojo Rubio Margin
## 0.068 0.180 0.045 0.080 0.372
## 0.029 0.076 0.019 0.034 0.157
## 0.020 0.052 0.013 0.023 0.108
## 0.066 0.175 0.044 0.078 0.363
## Margin 0.182 0.483 0.120 0.215 1.000
##
## $dev.frequencies
##
## Negro Castano Rojo Rubio
## Oscuros 0.047 0.021 -0.001 -0.068
## Pardos -0.003 0.015 0.005 -0.017
## Verdes -0.011 -0.003 0.011 0.004
## Azules -0.032 -0.034 -0.015 0.081
##
## $dev.values
##
## Negro Castano Rojo Rubio
## Oscuros 27.865 12.716 -0.385 -40.196
## Pardos -1.966 9.071 2.846 -9.951
## Verdes -6.676 -1.919 6.324 2.270
## Azules -19.223 -19.868 -8.785 47.877
##
## $Entropies
## entropía relativa
## filas 1.82786 0.91393
## columnas 1.79823 0.89911
## conjunta 3.44765 0.86191
## esperada 3.62609 0.90652
##
## $Cond.entropies
## fila|columnas columna|filas
## 1.64942 1.61979
##
## $Information
## mutual_information G² statistics
## 0.17844 146.44358
##
## $Phi_stats
## estadísticas
## phi-cuadrado 0.23360
## chi-cuadrado 138.28984
## grados de libertad 9.00000
## C de Pearson 0.43516
## T de Tchuprow 0.27904
## phi de Cramer 0.27904