En esta sesión deseo que entiendas por qué es necesario ir más allá de la correlación (Pearson, Spearman, etc.) o las diferencias de valores centrales (t test, kruska wallis, etc.). Esta necesidad de ir más allá trae una técnica conocida como la regresión.

La sesión va en tres partes:

De Correlación a Regresión Regresión Lineal

#—– 1. De correlacion a regresion ——

Trabajemos con los datos del archivo en SPSS hsb.sav:

library(rio)
linkToData='https://github.com/PsicologiaPUCP/ArchivosDeDatos/raw/master/hsb_ok.sav'
hsb=import(linkToData)

Antes de correr cualquier análisis estadístico, debes revisar como el tipo de datos que tu archivo trae y :

str(hsb)
## 'data.frame':    600 obs. of  15 variables:
##  $ ID    : chr  "    1" "    2" "    3" "    4" ...
##   ..- attr(*, "label")= chr "Student identification number"
##   ..- attr(*, "format.spss")= chr "A5"
##  $ SEX   : num  2 1 2 2 2 1 1 2 1 2 ...
##   ..- attr(*, "label")= chr "Sex"
##   ..- attr(*, "format.spss")= chr "F5.0"
##   ..- attr(*, "labels")= Named num [1:2] 1 2
##   .. ..- attr(*, "names")= chr [1:2] "Male" "Female"
##  $ RACE  : num  2 2 2 2 2 2 2 2 2 2 ...
##   ..- attr(*, "label")= chr "Race or ethnicity"
##   ..- attr(*, "format.spss")= chr "F5.0"
##   ..- attr(*, "labels")= Named num [1:4] 1 2 3 4
##   .. ..- attr(*, "names")= chr [1:4] "Hispanic" "Asian" "Black" "White"
##  $ SES   : num  1 1 1 2 2 2 1 1 2 1 ...
##   ..- attr(*, "label")= chr "Socioeconomic status"
##   ..- attr(*, "format.spss")= chr "F5.0"
##   ..- attr(*, "labels")= Named num [1:3] 1 2 3
##   .. ..- attr(*, "names")= chr [1:3] "Low" "Medium" "High"
##  $ SCTYP : num  1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "label")= chr "School type"
##   ..- attr(*, "format.spss")= chr "F5.0"
##   ..- attr(*, "labels")= Named num [1:2] 1 2
##   .. ..- attr(*, "names")= chr [1:2] "Public" "Private"
##  $ HSP   : num  3 2 2 3 3 2 1 1 1 1 ...
##   ..- attr(*, "label")= chr "High school program"
##   ..- attr(*, "format.spss")= chr "F5.0"
##   ..- attr(*, "labels")= Named num [1:3] 1 2 3
##   .. ..- attr(*, "names")= chr [1:3] "General" "Academic preparatory" "Vocational/Technical"
##  $ LOCUS : num  0.29 -0.42 0.71 0.06 0.22 0.46 0.44 0.68 0.06 0.05 ...
##   ..- attr(*, "label")= chr "Locus of control Standardized to mean 0 & SD of 1"
##   ..- attr(*, "format.spss")= chr "F5.2"
##  $ CONCPT: num  0.88 0.03 0.03 0.03 -0.28 0.03 -0.47 0.25 0.56 0.15 ...
##   ..- attr(*, "label")= chr "Self-concept Standardized to mean 0 & SD of 1"
##   ..- attr(*, "format.spss")= chr "F5.2"
##  $ MOT   : num  0.67 0.33 0.67 0 0 0 0.33 1 0.33 1 ...
##   ..- attr(*, "label")= chr "Motivation Average of 3 motivation items"
##   ..- attr(*, "format.spss")= chr "F5.2"
##  $ CAR   : num  10 2 9 15 1 11 10 9 9 11 ...
##   ..- attr(*, "label")= chr "Career choice"
##   ..- attr(*, "format.spss")= chr "F5.0"
##   ..- attr(*, "labels")= Named num [1:17] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..- attr(*, "names")= chr [1:17] "Clerical" "Craftsman" "Farmer" "Homemaker" ...
##  $ RDG   : num  33.6 46.9 41.6 38.9 36.3 49.5 62.7 44.2 46.9 44.2 ...
##   ..- attr(*, "label")= chr "Reading T-score standardized to mean 50 and SD 10"
##   ..- attr(*, "format.spss")= chr "F5.2"
##  $ WRTG  : num  43.7 35.9 59.3 41.1 48.9 46.3 64.5 51.5 41.1 49.5 ...
##   ..- attr(*, "label")= chr "Writing T-score standardized to mean 50 and SD 10"
##   ..- attr(*, "format.spss")= chr "F5.2"
##  $ MATH  : num  40.2 41.9 41.9 32.7 39.5 46.2 48 36.9 45.3 40.5 ...
##   ..- attr(*, "label")= chr "Math T-score standardized to mean 50 and SD 10"
##   ..- attr(*, "format.spss")= chr "F5.2"
##  $ SCI   : num  39 36.3 44.4 41.7 41.7 41.7 63.4 49.8 47.1 39 ...
##   ..- attr(*, "label")= chr "Science T-score standardized to mean 50 and SD 10"
##   ..- attr(*, "format.spss")= chr "F5.2"
##  $ CIV   : num  40.6 45.6 45.6 40.6 45.6 35.6 55.6 55.6 55.6 50.6 ...
##   ..- attr(*, "label")= chr "Science T-score standardized to mean 50 and SD 10"
##   ..- attr(*, "format.spss")= chr "F5.2"

Casi todo salió numerico, pero no sabremos qué ajustar si no leemos el manual metodológico o el diccionario de datos o la metadata disponible.

De ahi que depemos pre procesar:

# a nominal
categoricals=c("SEX","RACE","SES","SCTYP","HSP","CAR")

hsb[,categoricals]=factorize(hsb[,categoricals])


# a ordinal:
hsb$SES=ordered(hsb$SES,levels=c("Low","Medium","High" ))

La variable de interés: Asumamos que nuestra variable de interés es el desempeño en matemáticas; así, nuestra variable dependiente está representada por la variable MATH.

Apartir de ahí, consideremos que nos interesa saber la posible relación que pueda tener la variable que ha medido el desempeño en escritura; así, una variable independiente sería la representada por la variable WRTG. Hasta ahora sabemos que como son dos variables de tipo numérico debemos usar una correlación. La gráfica de correlación es esta:

library(ggplot2)

base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
scatter = base + geom_point()
scatter

Vemos que hay una aparente relación. Calculemos los indices de correlación:

f1=formula(~MATH + WRTG)

# camino parametrico
pearsonf1=cor.test(f1,data=hsb)[c('estimate','p.value')]

# camino no parametrico
spearmanf1=cor.test(f1,data=hsb,method='spearman')[c('estimate','p.value')]
## Warning in cor.test.default(x = structure(c(40.2, 41.9, 41.9, 32.7, 39.5, :
## Cannot compute exact p-value with ties

Asumiendo un camino paramétrico, podemos pedir el coeficiente de Pearson, el cuál al ser calculado obtenemos 0.6326664 (con p-value= 0).

Si hubieramos seguido una ruta no paramétrica, informaríamos el coeficiente de Spearman:0.6415126 (con p-value= 0).

Ahora, consideremos que nos interesa saber a la vez la posible relación que pueda tener la variable que ha medido el desempeño en ciencias; así otra variable independiente sería la representada por la variable SCI. Como es otra variable numérica no podemos calcular la correlación de tres variables, pero podemos tratar de verlo visualmente:

library(scatterplot3d)
scatterplot3d(hsb[,c('SCI','WRTG','MATH')])

base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
base + geom_point(aes(color = SCI))

Calulemos las correlaciones:

f2=formula(~MATH+SCI)

# camino parametrico
pearsonf2=cor.test(f2,data=hsb)[c('estimate','p.value')]

# camino no parametrico
spearmanf2=cor.test(f2,data=hsb, method='spearman')[c('estimate','p.value')]
## Warning in cor.test.default(x = structure(c(40.2, 41.9, 41.9, 32.7, 39.5, :
## Cannot compute exact p-value with ties

Podríamos calcular la correlación de SCI con MATH, obteniendo el Pearson (0.6495261, p-value= 0) y el Spearman (0.6551515,p-value= 0).

Visualmente vemos relación, pero no tenemos un coeficiente para medir ello.

Finalmente, ¿si quisiéramos ver si el sexo tiene algun rol en todo esto? Como ésta es una variable categórica y dicotómica, lo primero que puede venir a la mente es esta gráfica:

base=ggplot(data=hsb, aes(x=SEX, y=MATH))
base + geom_boxplot(notch = T)

Los boxplots tienen un notch flanqueando a la mediana, para sugerir igualdad de medianas si éstos se intersectan; de ahi que parece no haber diferencia sustantiva entre hombres y mujeres en cuanto a su desempeño en MATH.

Una alternativa al boxplot seria las barras de error:

library(ggpubr)
ggerrorplot(data=hsb, x = "SEX", y = "MATH")

En este último caso, si las lineas (denotado por las barras de error de la media) se intersectan, sugeriria que los valores medios (en este caso la media) podrian ser iguales.

Verificar si hay o no igualdad entre distribuciones depende si las variables se distribuyen o no de manera normal:

library(ggplot2)
ggplot(hsb,aes(x=MATH)) + 
  geom_histogram(aes(y = ..density..),bins = 20, fill='green') +
  stat_function(fun = dnorm, colour = "red",
                      args = list(mean = mean(hsb$MATH, na.rm = TRUE),
                                 sd = sd(hsb$MATH, na.rm = TRUE))) + 
  facet_grid(~SEX) + 
  coord_flip()

Nota que los histogramas de la data real tienen encima la curva normal que idealmente tendría esa data. La lejanía entre ellos, sugeriría no normalidad.

Se suele usar un qqplot para explorar la presencia/ausencia de normalidad:

# se sugiere normalidad si los puntos no se alejan de la diagonal.

library(ggpubr)
ggqqplot(data=hsb,x="MATH") + facet_grid(. ~ SEX)

Como ello no es fácil de discernir visualmente, tenemos por costumbre calcular algun coeficiente, como el Shapiro-Wilk:

library(knitr)
library(magrittr)
library(kableExtra)
f4=formula(MATH~SEX)


tablag= aggregate(f4, hsb,
          FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})

# para que se vea mejor:


shapiroTest=as.data.frame(tablag[,2])
names(shapiroTest)=c("W","Prob")


kable(cbind(tablag[1],shapiroTest))%>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F, position = "left")
SEX W Prob
Male 0.9837903 0.0034565
Female 0.9790040 0.0001031

Esto nos sugiere un camino no paramétrico para ver la diferencia de valores medios, por lo que deberiamos usar la prueba de Mann-Whitney en vez de la prueba t para testaer la relación entre ambas.

tf4=t.test(f4,data=hsb)[c('estimate','p.value')]
wilcoxf4=wilcox.test(f4,data=hsb)['p.value']

La prueba no paramétrica no rechazaría la igualdad de valores medios (Mann-Whitney con p valor = 0.3085543).

base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
base + geom_point(aes(color = SEX))

base + geom_point(aes(size = SCI, color=SEX,shape=SEX)) + scale_shape_manual(values=c(3, 2))

Otra alternativa puede ser:

base + geom_point(aes(color = SCI)) + facet_grid(~SEX)

Y claro:

paleta <- c("coral1","cyan" )
colors <- paleta[as.numeric(hsb$SEX)]
scatterplot3d(hsb[,c('SCI','WRTG','MATH')],color=colors)

En todos los gráficos vemos que los hombres y las mujeres están distribuidos por todo el gráfico, lo cual nos sugiere que no hay diferencias aun en dimensiones mayores a dos. Sin embargo, no tenemos una medida de cuanto cada uno afecta a nuestra dependiente.

De ahi que necesitamos la regresión.

#———- 2. Regresión Lineal————-

La regresión es una técnica donde hay que definir una variable dependiente y una o más independientes. Las independientes pueden tener rol predictor, dependiendo del diseño de investigación, aunque por defecto tiene un rol explicativo.

La regresión sí quiere informar cuánto una variable (independiente) puede explicar la variación de otra (dependiente), de ahí que es una técnica para probar hipótesis direccionales o asimétricas (las correlaciones tiene hipótesis simétricas).

La regresión busca proponer un modelo, es decir una ecuación, que recoja como una variable explicaría a otra.Para nuestro caso la variable dependiente es MATH, y tendriamos hasta ahora tres modelos:

modelo1=formula(MATH~WRTG)
modelo2=formula(MATH ~ WRTG + SCI)
modelo3= formula(MATH ~ WRTG + SCI + SEX)

Por ejemplo, para la hipótesis ‘el nivel de desempeño en escritura afecta el desempeño en matemáticas’, la regresión arrojaría este resultado:

library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
reg1=lm(modelo1,data=hsb)
stargazer(reg1,type = "html",intercept.bottom = FALSE)
## 
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>MATH</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Constant</td><td>19.769<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(1.633)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">WRTG</td><td>0.612<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.031)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>600</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.400</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.399</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>7.297 (df = 598)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>399.110<sup>***</sup> (df = 1; 598)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Aquí ya sabemos algo interesante, primero que WRTG tiene efecto, pues es significativo (indicado por los dos asteriscos); segundo, que ese efecto es directo, pues el coeficiente calculado es positivo; y tercero que la magnitud de ese efecto es 0.612, lo que indica cuanto aumenta MATH en promedio cuando WRTG se incremente en una unidad.

Esto es información suficiente para representar esa relación con una ecuación. Como la ecuación sólo tiene una variable independiente podemos producir una recta sobre el gráfico de correlación:

ggplot(hsb, aes(x=WRTG, y=MATH)) + 
  geom_point()+
  geom_smooth(method=lm)
## `geom_smooth()` using formula 'y ~ x'

Esa recta podemos representarla así:

El Y verdadero es MATH, pero la regresión produce un MATH^ estimado, de ahi la presencia del ϵ. Justamente el R cuadrado ajustado (0.4002668) nos brinda un porcentaje (multiplicalo por 100) que da una pista de nuestra cercanía a una situación perfecta (cuando vale 1).

Y sí queremos ver el efecto de SCI?

reg2=lm(modelo2,data=hsb)
stargazer(reg2,type = "html",intercept.bottom = FALSE)
## 
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>MATH</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Constant</td><td>10.629<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(1.630)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">WRTG</td><td>0.377<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.033)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">SCI</td><td>0.415<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.033)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>600</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.524</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.523</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>6.505 (df = 597)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>328.846<sup>***</sup> (df = 2; 597)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

En este caso, la regresión tendrá una formula con dos variables explicando la dependiente, así que en vez de producir una línea buscará producir un plano:

library(scatterplot3d)
G  <- scatterplot3d(hsb[,c('SCI','WRTG','MATH')])
G$plane3d(reg2, draw_polygon = TRUE, draw_lines = FALSE)

Este plano podemos representarlo así:

Nuevamente, el Y verdadero es MATH, pero la regresión produce un MATH^ estimado en forma de plano. De igual manera el R cuadrado ajustado (0.5241861) nos da una pista de nuestra lejanía a una situación perfecta.

Es clave darse cuenta de otro detalle, que el coeficiente de WRTG ha variado en la fórmula ahora que está presente SCI ¿Por qué sucede esto? Veamoslo así: en el primer caso, WRTG y ϵ buscaban representar la variabilidad en MATH, y ahora, en el segundo caso, viene SCI para mejorar esa explicación; es así que el peso de la explicación ahora se recalcula y el coeficiente de WRTG deja de explicar lo que le corresponde a SCI, y ϵ también le entrega algo a SCI.

Como ϵ no tiene coeficiente, representamos su variación usando el error típico de los residuos o residual standard error (RSE). Nótese que éste ha variado de un modelo ha otro, ahora tenemos un RSE menor. Aquí vale la pena preguntarse si esta disminución del error es significativa, obteniendo:

tanova=anova(reg1,reg2)
stargazer(tanova,type = 'html',summary = F,title = "Table de Análisis de Varianza")
## 
## <table style="text-align:center"><caption><strong>Table de Análisis de Varianza</strong></caption>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>Res.Df</td><td>RSS</td><td>Df</td><td>Sum of Sq</td><td>F</td><td>Pr(> F)</td></tr>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">1</td><td>598</td><td>31,842.070</td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">2</td><td>597</td><td>25,262.730</td><td>1</td><td>6,579.335</td><td>155.481</td><td>0</td></tr>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr></table>

La comparación de modelos usando la tabla de análisis de varianza (anova) propone como hipótesis nula que los modelos no difieren (no se ha reducido el error al pasar de un modelo a otro). Como la comparación es significativa (vea el P), rechazamos igualdad de modelos: el modelo 2 sí reduce el error al incluir una variable más.

Finalmente, veamos el rol de sexo:

reg3=lm(modelo3,data=hsb)
stargazer(reg3,type = "html",intercept.bottom = FALSE)
## 
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>MATH</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Constant</td><td>11.306<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(1.630)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">WRTG</td><td>0.424<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.036)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">SCI</td><td>0.375<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.035)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">SEXFemale</td><td>-1.922<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.582)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>600</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.533</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.530</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>6.452 (df = 596)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>226.510<sup>***</sup> (df = 3; 596)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Aunque no podemos graficar cuatro coordenadas, podemos usar elementos visuales:

library(scatterplot3d)
colors <- paleta[as.numeric(hsb$SEX)]
G  <- scatterplot3d(hsb[,c('SCI','WRTG','MATH')],color=colors)
G$plane3d(reg2, draw_polygon = TRUE, draw_lines = FALSE)

Nuestra nueva ecuación sería:

Nuevamente podemos ver si añadir SEXO en este modelo representa una mejora al anterior:

tanova=anova(reg1,reg2,reg3)
stargazer(tanova,type = 'html',summary = F,title = "Table de Análisis de Varianza")
## 
## <table style="text-align:center"><caption><strong>Table de Análisis de Varianza</strong></caption>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td>Res.Df</td><td>RSS</td><td>Df</td><td>Sum of Sq</td><td>F</td><td>Pr(> F)</td></tr>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">1</td><td>598</td><td>31,842.070</td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">2</td><td>597</td><td>25,262.730</td><td>1</td><td>6,579.335</td><td>158.063</td><td>0</td></tr>
## <tr><td style="text-align:left">3</td><td>596</td><td>24,808.420</td><td>1</td><td>454.308</td><td>10.914</td><td>0.001</td></tr>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr></table>

Finalmente, podemos resumir todo en esta tabla:

library(stargazer)
stargazer(reg1,reg2,reg3, type = "html", title = "Modelos planteadas",digits = 2, single.row = F,no.space = F,intercept.bottom = FALSE,
          dep.var.caption="Variable dependiente:",
          dep.var.labels="Desempeño en Matemáticas",
          covariate.labels=c("Constante","Desempeño en Escritura","Desempeño en Ciencias","SEXO (mujer)"),
          keep.stat = c("n","adj.rsq","ser"),df = F,
          notes.label = "Notas:")
## 
## <table style="text-align:center"><caption><strong>Modelos planteadas</strong></caption>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="3">Variable dependiente:</td></tr>
## <tr><td></td><td colspan="3" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="3">Desempeño en Matemáticas</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Constante</td><td>19.77<sup>***</sup></td><td>10.63<sup>***</sup></td><td>11.31<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(1.63)</td><td>(1.63)</td><td>(1.63)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Desempeño en Escritura</td><td>0.61<sup>***</sup></td><td>0.38<sup>***</sup></td><td>0.42<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.03)</td><td>(0.03)</td><td>(0.04)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Desempeño en Ciencias</td><td></td><td>0.42<sup>***</sup></td><td>0.37<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(0.03)</td><td>(0.04)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">SEXO (mujer)</td><td></td><td></td><td>-1.92<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td>(0.58)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>600</td><td>600</td><td>600</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.40</td><td>0.52</td><td>0.53</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>7.30</td><td>6.51</td><td>6.45</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Notas:</td><td colspan="3" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Gráficamente:

library(ggplot2)
library(sjPlot)
## Registered S3 methods overwritten by 'parameters':
##   method                           from      
##   as.double.parameters_kurtosis    datawizard
##   as.double.parameters_skewness    datawizard
##   as.double.parameters_smoothness  datawizard
##   as.numeric.parameters_kurtosis   datawizard
##   as.numeric.parameters_skewness   datawizard
##   as.numeric.parameters_smoothness datawizard
##   print.parameters_distribution    datawizard
##   print.parameters_kurtosis        datawizard
##   print.parameters_skewness        datawizard
##   summary.parameters_kurtosis      datawizard
##   summary.parameters_skewness      datawizard
plot_models(reg1,reg2,reg3,vline.color = "grey",m.labels=c("Modelo 1","Modelo 2","Modelo 3"))

Como vemos, el propósito del último gráfico es mostrar que dado que ninguna de las variables (y sus intervalos de confianza) tocan el valor cero (que significa que la variable no tiene efecto en la dependiente).