Modelos estádisticos tradicionales como regresiones

Objetivo del manual

• Entender los abordajes estadísticos tradicionales como regresiones

 

Aquí veremos las pruebas estadísticas más comunes y mostraremos cómo pueden representarse en el formato de regresión lineal. Esta sección se basa en este artículo. Consúltelo para obtener una descripción más detallada de las alternativas no paramétricas.

---

1 Prueba t de una muestra

La prueba evalúa si la media de una variable continua es diferente de 0. Es equivalente a una regresión lineal sin predictor, que prueba si el intercepto es diferente de 0:

\(y = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\)

 

Código

# number of observations
n <- 50

# set seed
set.seed(123)

# create variable with mean 1
y <- rnorm(n = n, mean = 1)

# run t test
(t <- t.test(y))

    One Sample t-test

data:  y
t = 7.9, df = 49, p-value = 2.7e-10
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 0.77127 1.29753
sample estimates:
mean of x 
   1.0344 

Código

# run equivalent linear regression
(lm_t <- summary(lm(y ~ 1)))

Call:
lm(formula = y ~ 1)

Residuals:
   Min     1Q Median     3Q    Max 
-2.001 -0.594 -0.107  0.664  2.135 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    1.034      0.131     7.9  2.7e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.926 on 49 degrees of freedom

 

Pongamos juntos los resultados de ambas pruebas para compararlos más de cerca:

	model	estimate	p.value	t
mean of x	t-test	1.0344	0	7.9
	lm	1.0344	0	7.9

 

Tenga en cuenta que, como no hay predictores en el modelo, utilizamos un ‘1’ en el lugar de los predictores en la fórmula del modelo (‘y ~ 1’).

---

2 Prueba de t pareada

Una prueba t pareada evalúa si la diferencia de promedios entre dos variables numéricas es 0

\(y1 - y2 = \beta_0 \qquad \mathcal{H}_0: \beta_0 = 0\)

 

El modelo lineal correspondiente es el mismo que el de la prueba t de una muestra, pero la variable de entrada es la diferencia entre las dos variables (y1 - y2):

Código

set.seed(123)

# sample size
n <- 50

# variable with mean 1
y1 <- rnorm(n = n, mean = 1)

# variable with mean 3
y2 <- rnorm(n = n, mean = 1.4)

# run paired t test
paired_t <- t.test(y1, y2, paired = TRUE)

paired_t

    Paired t-test

data:  y1 and y2
t = -2.75, df = 49, p-value = 0.0084
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.88659 -0.13742
sample estimates:
mean of the differences 
                 -0.512 

Código

# difference between the 2 variables
diff_y <- y1 - y2

# run model
lm_paired_t <- summary(lm(formula = diff_y ~ 1))

lm_paired_t

Call:
lm(formula = diff_y ~ 1)

Residuals:
   Min     1Q Median     3Q    Max 
-2.600 -0.701  0.112  1.035  2.909 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   -0.512      0.186   -2.75   0.0084 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.32 on 49 degrees of freedom

	model	estimate	p.value	t
mean of the differences	paired t-test	-0.512	0.00839	-2.7468
	lm	-0.512	0.00839	-2.7468

 

---

3 Prueba t de dos promedios

También llamada prueba t independiente. Evalúa si los promedios de las dos variables son diferentes. La hipótesis nula es que la diferencia es 0:

\(y = \beta_0 + \beta_1 * x1 \qquad \mathcal{H}_0: \beta_1 = 0\)

Código

# set seed
set.seed(123)

# number of observations
n <- 50
b0 <- -4
b1 <- 3
error <- rnorm(n = n, mean = 0, sd = 1)

# random variables
x1_num <- rbinom(n = n, size = 1, prob = 0.5)
y <- b0 + b1 * x1_num + error

x1 <- factor(x1_num, labels = c("a", "b"))

# create data frame
xy_data_cat <- data.frame(x1, x1_num, y)

# run paired t test
indep_t <- t.test(xy_data_cat$y[xy_data_cat$x1 == "a"], xy_data_cat$y[xy_data_cat$x1 == "b"])

indep_t

    Welch Two Sample t-test

data:  xy_data_cat$y[xy_data_cat$x1 == "a"] and xy_data_cat$y[xy_data_cat$x1 == "b"]
t = -12.3, df = 47.9, p-value <2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.7518 -2.6957
sample estimates:
mean of x mean of y 
 -4.07748  -0.85371 

Código

# run regression model
lm_indep_t <- summary(lm(formula = y ~ x1, data = xy_data_cat))

lm_indep_t

Call:
lm(formula = y ~ x1, data = xy_data_cat)

Residuals:
   Min     1Q Median     3Q    Max 
-2.113 -0.594 -0.136  0.552  2.246 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -4.077      0.186   -22.0   <2e-16 ***
x1b            3.224      0.263    12.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.928 on 48 degrees of freedom
Multiple R-squared:  0.758, Adjusted R-squared:  0.753 
F-statistic:  151 on 1 and 48 DF,  p-value: <2e-16

	model	estimate	p.value	t
mean of y	2 means t-test	3.2238	0	-12.276
	lm	3.2238	0	12.276

 

---

4 Tres o más promedios: ANDEVA unidireccional

Muy similar a la prueba t de dos promedios, pero con tres o más niveles en la variable categórica:

\(\hat{Y} \sim \beta_{o} + \beta_{1} * x_{1} + \beta_{2} * x_{2} + ... \qquad \mathcal{H}_0: y = \beta_0\)

 

Se trata de un modelo de regresión múltiple. La variable categórica se codifica con dummy, por lo que se divide en predictores indicadores (\(x_i\) es \(x=0\) o \(x=1\)).

Un conjunto de datos con una variable categórica de 3 niveles puede simularse así:

Código

# set seed
set.seed(123)

# number of observations
n <- 50
b0 <- -4
b1 <- 3
error <- rnorm(n = n, mean = 0, sd = 1)

# random variables
x1_num <- rbinom(n = n, size = 2, prob = c(0.33, 0.33))
y <- b0 + b1 * x1_num + error

x1 <- factor(x1_num, labels = c("a", "b", "c"))

# create data frame
xy_data_cat2 <- data.frame(x1, y)

head(xy_data_cat2)

x1	y
b	-1.56048
a	-4.23018
b	0.55871
c	2.07051
b	-0.87071
b	0.71506

Código

# ANOVA function
anv_1w <- anova(aov(formula = y1 ~ x1, data = xy_data_cat2))

anv_1w

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
x1	2	4.3728	2.18638	2.7307	0.07552
Residuals	47	37.6318	0.80068	NA	NA

Código

# linear regression
lm_anv_1w <- summary(lm(formula = y1 ~ x1, data = xy_data_cat2))

lm_anv_1w

Call:
lm(formula = y1 ~ x1, data = xy_data_cat2)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6721 -0.5195 -0.0289  0.4839  2.3266 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.842      0.187    4.51  4.2e-05 ***
x1b            0.565      0.277    2.04    0.047 *  
x1c           -0.141      0.367   -0.38    0.703    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.895 on 47 degrees of freedom
Multiple R-squared:  0.104, Adjusted R-squared:  0.066 
F-statistic: 2.73 on 2 and 47 DF,  p-value: 0.0755

	model	F.statistic	p.value
	1 way anova	2.7307	0.07552
value	lm	2.7307	0.07552

 

---

5 Chi-cuadrado para tablas de contingencia

Esta prueba busca la asociación entre dos variables categóricas basándose en la co-ocurrencia de las categorías, medidas como recuentos. Esto es lo que es una tabla de contingencia: recuentos de todas las combinaciones de categorías entre dos variables. Por ejemplo, las muertes por edad de los pasajeros del Titanic:

Código

# Contingency table
# modified from https://lindeloev.github.io/tests-as-linear/
mood_data <- matrix(
  data = c(100, 70, 30, 32, 110, 120), nrow = 2,
  dimnames = list(handedness = c("left_handed", "right_handed"), mood = c("happy", "meh", "sad"))
)

mood_data

              mood
handedness     happy meh sad
  left_handed    100  30 110
  right_handed    70  32 120

 

Aunque la codificación es un poco más elaborada, puede reducirse a un modelo ANDEVA log-lineal de dos vías:

\(log(y_i) = log(N) + log(\alpha_i) + log(\beta_j) + log(\alpha_i\beta_j)\)

 

Los coeficientes de regresión son \(A_i\) y \(B_i\). \(\alpha_i\) y \(\beta_j\) son las proporciones de recuentos para los dos niveles de la variable categórica binaria (uno para cada uno de los niveles de la segunda variable categórica). Afortunadamente, existe una implementación de este modelo en el paquete R ‘MASS’ que simplifica la codificación:

Código

# Built-in chi-square. It requires matrix format.
chi_mod <- chisq.test(mood_data)

chi_mod

    Pearson's Chi-squared test

data:  mood_data
X-squared = 5.1, df = 2, p-value = 0.078

Código

# convert to long
mood_data_long <- as.data.frame(as.table(mood_data))

# log linear model from MASS package
log_lm_mod <- MASS::loglm(Freq ~ handedness + mood, data = mood_data_long)

log_lm_mod

Call:
MASS::loglm(formula = Freq ~ handedness + mood, data = mood_data_long)

Statistics:
                    X^2 df P(> X^2)
Likelihood Ratio 5.1199  2 0.077308
Pearson          5.0999  2 0.078087

	model	X.squared	p.value
X-squared	chi-square	5.0999	0.07809
	log lm	5.0999	0.07809

 

6 Alternativas “no paramétricas” como modelos de regresión lineal

Sólo presentaremos un ejemplo de prueba “no paramétrica”. No obstante, este único ejemplo debería bastar para demostrar la lógica que subyace a los procedimientos “no paramétricos”. La mayoría de las pruebas “no paramétricas” se basan en un sencillo truco de transformación de datos: la transformación de rango de las respuestas. Las transformaciones de rango convierten los valores en un índice según su posición cuando se ordenan por magnitud:

Código

# simulate variable
set.seed(123)
y <- round(rnorm(n = 5), 2)

y

[1] -0.56 -0.23  1.56  0.07  0.13

Código

# convert to ranks
rank_y <- rank(y)

rank_y

[1] 1 2 5 3 4

 

6.1 Prueba de Wilcoxon para dos promedios

También conocida como prueba de Mann-Whitney, es la alternativa no paramétrica a la prueba t de dos promedios:

Código

# set seed
set.seed(123)

# number of observations
n <- 50
y <- rnorm(n = n, mean = 0.2, sd = 1)

# create data frame
y_data <- data.frame(y)


y

 [1] -0.3604756 -0.0301775  1.7587083  0.2705084  0.3292877  1.9150650
 [7]  0.6609162 -1.0650612 -0.4868529 -0.2456620  1.4240818  0.5598138
[13]  0.6007715  0.3106827 -0.3558411  1.9869131  0.6978505 -1.7666172
[19]  0.9013559 -0.2727914 -0.8678237 -0.0179749 -0.8260044 -0.5288912
[25] -0.4250393 -1.4866933  1.0377870  0.3533731 -0.9381369  1.4538149
[31]  0.6264642 -0.0950715  1.0951257  1.0781335  1.0215811  0.8886403
[37]  0.7539177  0.1380883 -0.1059627 -0.1804710 -0.4947070 -0.0079173
[43] -1.0653964  2.3689560  1.4079620 -0.9231086 -0.2028848 -0.2666554
[49]  0.9799651  0.1166309

Código

# Wilcoxon / Mann-Whitney U
wilcx_mod <- wilcox.test(y_data$y)

wilcx_mod

    Wilcoxon signed rank test with continuity correction

data:  y_data$y
V = 797, p-value = 0.12
alternative hypothesis: true location is not equal to 0

Código

signed_rank <- function(x) sign(x) * rank(abs(x))

# As linear model with our dummy-coded group_y2:
wicx_lm_mod <- lm(signed_rank(y) ~ 1, data = y_data) # compare to

summary(wicx_lm_mod)

Call:
lm(formula = signed_rank(y) ~ 1, data = y_data)

Residuals:
   Min     1Q Median     3Q    Max 
-53.38 -24.13   0.12  25.37  43.62 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)     6.38       4.09    1.56     0.12

Residual standard error: 28.9 on 49 degrees of freedom

model	p.value
1 mean wilcoxon	0.12482
wilcoxon_lm	0.12480

 

---

7 Referencias

• Common statistical tests are linear models

 

---

 

Información de la sesión

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=es_CR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=es_CR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=es_CR.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] viridis_0.6.2     viridisLite_0.4.2 ggplot2_3.4.4     knitr_1.45       

loaded via a namespace (and not attached):
 [1] pillar_1.9.0      compiler_4.1.2    tools_4.1.2       digest_0.6.33    
 [5] jsonlite_1.8.7    evaluate_0.23     lifecycle_1.0.4   tibble_3.2.1     
 [9] gtable_0.3.4      pkgconfig_2.0.3   rlang_1.1.2       cli_3.6.1        
[13] DBI_1.1.3         rstudioapi_0.13   yaml_2.3.7        xfun_0.41        
[17] fastmap_1.1.1     gridExtra_2.3     withr_2.5.2       dplyr_1.0.8      
[21] generics_0.1.2    vctrs_0.6.4       htmlwidgets_1.5.4 grid_4.1.2       
[25] tidyselect_1.1.2  glue_1.6.2        R6_2.5.1          fansi_1.0.5      
[29] rmarkdown_2.25    purrr_0.3.4       magrittr_2.0.3    MASS_7.3-55      
[33] scales_1.2.1      htmltools_0.5.7   assertthat_0.2.1  colorspace_2.1-0 
[37] utf8_1.2.4        munsell_0.5.0