Distribuciones y estadísticas

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

Datos cualitativos

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