library(nnet)
library(pander)
## Warning: package 'pander' was built under R version 4.3.3

logit con factores inusuales

El siguiente ejemplo, realizado con una muestra intencional, tiene como factores significativos los grupos y no los grados; por lo que es cuestionable la existencia de factores inusuales. Además, el hecho de mezclar grupos de grados semejantes, hace que desaparezca la significatividad, lo cual se aleja de la teoría de efectos causales, que involucran un factor como la edad, por ejemplo, y no el hecho de hallar significatividad dentro de los grupos elegidos.

Este ejemplo, muestra que el modelo de regresión logística falla debido a que la muestra no es aleatoria.

p05: frecuencia semanal de estudio en casa (Nunca, casi nunca, a veces, casi siempre, siempre)

p02: frecuencia semanal de lectura de libros (Nunca, casi nunca, a veces, casi siempre, siempre)

p06: el estudiante está en 7B.

p07: el estudiante está en 8A.

P08: el estudiante está en 8B.

Se hace un modelo de regresión lineal múltiple para visualizar variables relacionadas con la respuesta. Se desea predecir el nivel de estudio en casa.

data0=read.delim('clipboard')
data0
##    id grado p01 p02 p03 p04 p05 p06 p07 p08
## 1   1    7B   0   2   2   3   2   1   0   0
## 2   2    7B   0   1   2   3   4   1   0   0
## 3   3    7B   0   2   2   3   2   1   0   0
## 4   4    7B   2   3   1   1   2   1   0   0
## 5   5    7B   0   4   2   4   2   1   0   0
## 6   6    8A   0   3   4   4   1   0   1   0
## 7   7    8B   1   2   4   4   4   0   0   1
## 8   8    8B   0   4   3   4   1   0   0   1
## 9   9    6A   0   4   1   1   4   0   0   0
## 10 10    9B   0   2   2   4   4   0   0   0
## 11 11    9B   0   2   2   4   4   0   0   0
## 12 12    8A   0   1   2   2   0   0   1   0
## 13 13    8B   0   2   2   2   2   0   0   1
## 14 14    8A   2   2   1   3   2   0   1   0
## 15 15    8A   0   2   0   4   2   0   1   0
## 16 16    8B   0   2   4   2   2   0   0   1
## 17 17    8B   0   0   4   1   4   0   0   1
## 18 18    8B   0   4   4   2   2   0   0   1
## 19 19    8A   0   2   1   4   2   0   1   0
## 20 20    8A   0   1   2   3   2   0   1   0
## 21 21    8A   0   0   2   2   3   0   1   0
## 22 22    8A   0   2   0   1   4   0   1   0
lre=lm(data0$p05~data0$p02+data0$p06+data0$p07+data0$p08)
summary(lre)
## 
## Call:
## lm(formula = data0$p05 ~ data0$p02 + data0$p06 + data0$p07 + 
##     data0$p08)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.26919 -0.53116  0.00997  0.29736  2.16152 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.1486     0.7631   6.746 3.42e-06 ***
## data0$p02    -0.4307     0.1940  -2.220  0.04032 *  
## data0$p06    -1.7149     0.7114  -2.410  0.02754 *  
## data0$p07    -2.4487     0.6881  -3.558  0.00242 ** 
## data0$p08    -1.6436     0.6901  -2.382  0.02918 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9716 on 17 degrees of freedom
## Multiple R-squared:  0.456,  Adjusted R-squared:  0.328 
## F-statistic: 3.562 on 4 and 17 DF,  p-value: 0.02759
head(probability.table<-fitted(lre))
##        1        2        3        4        5        6 
## 2.572283 3.002991 2.572283 2.141575 1.710867 1.407777

A continuación se obtienen las probabilidades de predicción de cada grupo.

data0=read.delim('clipboard')
beta=cbind(c(-133.00976,59.39987, -29.51153,15.6998408,24.083102),c(-33.71131,20.98523,85.95106,12.0150243,78.433796),c(17.64231,-44.42444,14.24601,-0.9522501,-5.685661),c(109.13267,18.98371,-54.42722,-128.7921263,-61.063578))
data1=rbind(c(1,    2,  1,  0,  0), c(1,    1,  1,  0,  0 ), c(1,   2,  1,  0,  0), c(1,    3,  1,  0,  0), c(1,    4,  1,  0,  0),c(1, 3,  0,  1,  0), c(1,    2,  0,  0,  1), c(1,    4,  0,  0,  1), c(1,    4,  0,  0,  0), c(1,    2,  0,  0,  0),c(1, 2,  0,  0,  0), c(1,    1,  0,  1,  0), c(1,    2,  0,  0,  1), c(1,    2,  0,  1,  0), c(1,    2,  0,  1,  0), c(1,    2,  0,  0,  1), c(1,    0,  0,  0,  1), c(1,    4,  0,  0,  1), c(1,    2,  0,  1,  0), c(1,    1,  0,  1,  0), c(1,    0,  0,  1,  0),c(1, 2,  0,  1,  0))
pm<-data1 %*% beta
em=exp(pm)
round(rbind(cbind(1/(1+apply(em,1,sum)),em/(1+apply(em,1,sum)))),8)
##             [,1]       [,2]       [,3]      [,4]       [,5]
##  [1,] 0.00000000 0.00000000 0.82307771 0.0000000 0.17692229
##  [2,] 0.00000000 0.00000000 0.38599468 0.0000000 0.61400532
##  [3,] 0.00000000 0.00000000 0.82307771 0.0000000 0.17692229
##  [4,] 0.00000000 0.00000000 0.97177346 0.0000000 0.02822654
##  [5,] 0.00000000 0.00000000 0.99609031 0.0000000 0.00390969
##  [6,] 0.00000000 1.00000000 0.00000000 0.0000000 0.00000000
##  [7,] 0.00000000 0.00000000 0.65845888 0.0000000 0.34154112
##  [8,] 0.00000000 0.49999728 0.49531140 0.0000000 0.00469131
##  [9,] 0.00000000 0.00000000 0.00000000 0.0000000 1.00000000
## [10,] 0.00000000 0.00000000 0.00000000 0.0000000 1.00000000
## [11,] 0.00000000 0.00000000 0.00000000 0.0000000 1.00000000
## [12,] 0.50002450 0.00000000 0.24557475 0.0000000 0.25440075
## [13,] 0.00000000 0.00000000 0.65845888 0.0000000 0.34154112
## [14,] 0.00000000 0.00000001 0.87720350 0.0000000 0.12279649
## [15,] 0.00000000 0.00000001 0.87720350 0.0000000 0.12279649
## [16,] 0.00000000 0.00000000 0.65845888 0.0000000 0.34154112
## [17,] 0.00000000 0.00000000 0.03400648 0.0000000 0.96599352
## [18,] 0.00000000 0.49999728 0.49531140 0.0000000 0.00469131
## [19,] 0.00000000 0.00000001 0.87720350 0.0000000 0.12279649
## [20,] 0.50002450 0.00000000 0.24557475 0.0000000 0.25440075
## [21,] 0.00000006 0.00000000 0.00000000 0.9999999 0.00000000
## [22,] 0.00000000 0.00000001 0.87720350 0.0000000 0.12279649

El porcentaje de clasificaciones correctas es $18/22*100%=81.82%.