En análisis de correlación y de regresión simple son técnicas para identificar y modelar o representar la relación que existe entre dos variables cuantitativas medidas en escala de intervalo o razón.

Se dice que dos variables están correlacionadas cuando se aprecia un cambio sistemático en sus respectivas puntuaciones.

El análisis de regresión simple busca identificar y representar una relación lineal entre dos variables de intervalo o razón. Una relación lineal puede representarse mediante una ecuación lineal o una recta de regresión:

\[\hat{Y}=b_0 + b_{1}X\]

Esta ecuación nos dice que el valor esperado de Y (\(\hat{Y}\)) será igual a una constante \(b_0\) más el puntaje de X multiplicado por un coeficiente \(b_1\). En tal sentido, por cada cambio de \(X+1\) se espera que \(\hat{Y}\) cambie en \(b_1\). Por otro lado, si \(X=0\), entonces \(\hat{Y}=b_0\)

Indicadores sociodemográficos para 207 países del mundo en 1998. Los datos se encuentran en el archivo "BD_Mundo.zip" que puede descargarse desde:

http://sites.google.com/a/pucp.pe/data_est/archivos

load("mundo98.rda")
## Cargamos además algunos paquetes que usaremos en esta clase:
library(ggplot2)
library(grid)
library(scales)

Vamos a trabajar con dos variables:

• educationFemale: Promedio de años de educación formal de las mujeres
• contraception: Porcentaje de mujeres que utilizan métodos anticonceptivos

myvars <- c("educationFemale", "contraception")
data <- na.omit(mundo98[myvars])

En el análisis de regresión usualmente debemos seguir los siguientes pasos:

1. Inspeccionar visualmente la relación entre las variables mediante un diagrama de dispersión.
2. Estimar la ecuación o recta de regresión.
3. Interpretar los coeficientes de regresión.
4. Calcular los coeficientes de correlación y de determinación
5. Evaluar la significancia de los coeficientes de regresión y correlación

p <- ggplot(data, aes(x=educationFemale, y=contraception)) + geom_point()
p

La recta o ecuación de regresión es una línea recta que pasa en medio de los puntos del diagrama de dispersión y que representa la relación entre ambas variables.

Para calcularla se emplea el criterio de los mínimos cuadrados

El método de los mínimos cuadrados implica determinar los valores de los coeficientes de regresión \(b_0\) y \(b_1\) que minimizan la suma de cuadrados de las desviaciones entre los valores observados de Y y sus respectivos valores estimados \(\hat{Y}\)

\[\sum{e^2}=\sum(Y_{i}-\hat{Y}_i)^2\]

\[e_{Peru}^2 = (Y_{i}-\hat{Y}_{i})^2=(64-53.56)^2=(10.44)^2=108.99\]

La recta que minimiza los errores cuadrados se puede calcular con las siguientes ecuaciones:

\[b_{1}=\frac{\sum(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sum(x_{i}-\bar{x})^2}\]

\[b_{0}=\bar{y}-b_{1}\bar{x}\]

En el R con la siguiente sintaxis se puede obtener los coeficientes de regresión \(b_0\) y \(b_1\)

lm(contraception~educationFemale, data=data)

## 
## Call:
## lm(formula = contraception ~ educationFemale, data = data)
## 
## Coefficients:
##     (Intercept)  educationFemale  
##          -9.461            5.341

Por lo tanto: \(\hat{Y}=-9.461 + 5.341(X)\)

De acuerdo con los resultados del cálculo de los coeficientes, nuestro modelo o recta de regresión simple nos dice que:

\[\hat{Y}=-9.461 + 5.341(X)\]

Eso significa que:

• Por cada incremento de 1 año en los años promedio de estudio de las mujeres en un país, se espera que el % de mujeres que usan anticonceptivos se incremente en 5.341%.
• Si el promedio de años de estudio de las mujeres en un país fuese = 0, se esperaría un % de uso de anticonceptivos de -9.461% (es un valor teórico)

Prueba para un caso: En el Perú, en 1998 en promedio las mujeres estudiaban 11.8 años. De acuerdo con el modelo calculado, dado ese valor de X se esperaría que el porcentaje de mujeres que usan anticonceptivos en el Perú fuese de 53.56%:

\[\hat{Y}=-9.461+(5.341)(11.8)=53.56\]

El coeficiente de correlación bivariable r de Pearson mide la fuerza y dirección de la relación entre dos variables cuantitativas en una escala que varía entre -1 y +1. Cuanto más alejado del 0 sea el valor del coeficiente, más fuerte será la relación. El signo nos indica si se trata de una relación directa o inversa.

El coeficiente r de Pearson se calcula utilizando la siguiente fórmula

\[r=\frac{\sum(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum(x_{i}-\bar{x})^2\sum(y_{i}-\bar{y})^2}}\]

En el R, el coeficiente de correlación se pide usando el siguiente comando:

cor(data$contraception, data$educationFemale)

## [1] 0.8192834

• El \(R^2\) mide la "bondad de ajuste" del modelo de regresión a los datos analizados.
• Nos indica qué tan bien la recta de regresión es capaz de predecir los valores de Y
• \(R^2\) varía en una escala de 0 a 1, cuanto mayor es su valor, mayor poder predictivo tiene el modelo de regresión.
• Puede interpretarse como la proporción de la varianza total de Y que es "explicada" por el modelo de regresión (similar a ANOVA)

En el caso de Perú

\[Y-\bar{Y}=(\hat{Y}-\bar{Y})+(Y-\hat{Y})\]

\[64-50.9=(53.56-50.9)+(64-53.56)\]

\[13.1=2.66+10.44\]

En el caso de Zuiza

\[71-50.9=(62.46-50.9)+(71-62.46)\]

\[20.1=11.56+8.54\]

Si aplicamos la misma lógica a todos los casos tenemos:

\[\sum(Y-\bar{Y})^2=\sum(\hat{Y}-\bar{Y})^2+\sum(Y-\hat{Y})^2\]

La suma total de cuadrados \(SS_y\) es igual a la suma de cuadrados de la regresión \(SS_x\) más la suma de los errores cuadrados \(SS_e\).

\(R^2\) se calcula:

\[R^2=\frac{\sum(\hat{Y}-\bar{Y})^2}{\sum(Y-\bar{Y})^2}\]

Por lo tanto representa la proporción de la varianza total de Y que es "captada" o explicada por el modelo de regresión (X).

El \(R^2\) forma parte del conjunto de estadísticos de resumen de un modelo de regresión.

modelo1 <- lm(contraception~educationFemale, data=data)
summary(modelo1)

## 
## Call:
## lm(formula = contraception ~ educationFemale, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.088  -7.949   1.531   8.560  27.628 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -9.4606     5.7131  -1.656    0.103    
## educationFemale   5.3412     0.4826  11.068 3.99e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.39 on 60 degrees of freedom
## Multiple R-squared:  0.6712, Adjusted R-squared:  0.6657 
## F-statistic: 122.5 on 1 and 60 DF,  p-value: 3.992e-16

Para cada coeficiente de regresión se puede calcular un error estándar del coeficiente.

\[\sigma_{b_0}=\sqrt{\frac{\sigma_{e}^2\sum{x^2}}{n\sum(x-\bar{x})^2}}\]

\[\sigma_{b_1}=\sqrt{\frac{\sigma_{e}^2}{\sum(x-\bar{x})^2}}\]

Donde:

\[\sigma_{e}^2=\frac{\sum{e^2}}{(n-1-k)}=\frac{\sum(y-\hat{y})^2}{(n-1-k)}\]

Siendo \((n-1-k)\) los grados de libertad del error o residuales, donde \(k\) es el número de variables independientes del modelo.

Con los errores estándar se puede calcular el intervalo de confianza de los coeficientes del modelo de regresión:

\[\beta_{1-\alpha}=b \pm t_{\alpha/2}\sigma_{b}\]

En el modelo calculado tenemos que:

\[\sigma_{b_0}=5.71\]

\[\sigma_{b_1}=0.48\]

Se puede usar para ello la siguiente función para un intervalo de confianza del 95%

confint(modelo1)

##                      2.5 %   97.5 %
## (Intercept)     -20.888497 1.967228
## educationFemale   4.375853 6.306491

ggplot(data, aes(x=educationFemale, y=contraception)) + geom_point() + 
  geom_smooth(method=lm)

Con el error estándar de los coeficientes se puede calcular una prueba de hipótesis (prueba de t) para determinar si el coeficiente de regresión es estadísticamente significativo. En este caso las hipótesis se formulan:

\[H0: b = 0\] \[H1: b \neq 0\]

El estadístico de la prueba es un valor de \(t\) que se calcula:

\[t=\frac{b}{\sigma_b}\]

En nuestro ejemplo, para:

• \(b_0\) tenemos que \(t=-1.656\) con una significancia ó p-value de 0.103
• \(b_1\) tenemos que \(t=11.068\) con una significancia ó p-value de 3.99e-16

Podemos decidir si rechazamos H0 comparando el nivel de significancia con el \(\alpha\) que hayamos fijado para la prueba.

Para evaluar la significancia del coeficiente de determinación utilizamos la tabla de ANOVA del modelo:

anova(modelo1)

## Analysis of Variance Table
## 
## Response: contraception
##                 Df Sum Sq Mean Sq F value    Pr(>F)    
## educationFemale  1  21972 21971.5   122.5 3.992e-16 ***
## Residuals       60  10762   179.4                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El error estándar del modelo (o error residual estándar) es el error típico que se cometería al tratar de predecir el valor de Y usando la ecuación de regresion:

\[\sigma=\sqrt{\frac{\sum(y-\hat{y})^2}{(n-1-k)}}\]

En nuestro ejemplo, el error estándar del modelo es 13.39. Si queremos predecir el % de mujeres que usan anticonceptivos a partir de una regresión sobre la base de los años de educación promedio que tienen las mujeres en un país, el error típico sería de \(\pm 13.39\)%.

## 
## Call:
## lm(formula = contraception ~ educationFemale, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.088  -7.949   1.531   8.560  27.628 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -9.4606     5.7131  -1.656    0.103    
## educationFemale   5.3412     0.4826  11.068 3.99e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.39 on 60 degrees of freedom
## Multiple R-squared:  0.6712, Adjusted R-squared:  0.6657 
## F-statistic: 122.5 on 1 and 60 DF,  p-value: 3.992e-16

Paso 1: Diagrama de disperión de la Tasa de fecundidad según uso de anticonceptivos

ggplot(mundo98, aes(x=contraception, y=tfr)) + geom_point()

modelo2 <- lm(tfr~contraception, data=mundo98)
summary(modelo2)

## 
## Call:
## lm(formula = tfr ~ contraception, data = mundo98)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3000 -0.4562  0.0617  0.6620  2.4620 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.978989   0.186531   32.05   <2e-16 ***
## contraception -0.056477   0.003778  -14.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.069 on 139 degrees of freedom
##   (66 observations deleted due to missingness)
## Multiple R-squared:  0.6165, Adjusted R-squared:  0.6138 
## F-statistic: 223.5 on 1 and 139 DF,  p-value: < 2.2e-16

misvars <- c("tfr", "contraception")
data2 <- na.omit(mundo98[misvars])
cor(data2)

##                     tfr contraception
## tfr            1.000000     -0.785205
## contraception -0.785205      1.000000

summary(modelo2)

## 
## Call:
## lm(formula = tfr ~ contraception, data = mundo98)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3000 -0.4562  0.0617  0.6620  2.4620 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    5.978989   0.186531   32.05   <2e-16 ***
## contraception -0.056477   0.003778  -14.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.069 on 139 degrees of freedom
##   (66 observations deleted due to missingness)
## Multiple R-squared:  0.6165, Adjusted R-squared:  0.6138 
## F-statistic: 223.5 on 1 and 139 DF,  p-value: < 2.2e-16

ggplot(data2, aes(x=contraception, y=tfr)) + geom_point() + 
  geom_smooth(method=lm)

misvars <- c("tfr", "contraception")
data2 <- na.omit(mundo98[misvars])
cor.test(data2$tfr, data2$contraception)

## 
##  Pearson's product-moment correlation
## 
## data:  data2$tfr and data2$contraception
## t = -14.9498, df = 139, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8413116 -0.7123600
## sample estimates:
##       cor 
## -0.785205

Podemos evaluar rápidamente cómo se presentan las correlaciones entre un conjunto de variables utilizando una matriz de correlaciones:

misvars <- c("tfr", "contraception", "economicActivityFemale")
data3 <- na.omit(mundo98[misvars])
cor(data3)

##                               tfr contraception economicActivityFemale
## tfr                     1.0000000    -0.7490240             -0.2567071
## contraception          -0.7490240     1.0000000              0.1455538
## economicActivityFemale -0.2567071     0.1455538              1.0000000

library(Hmisc)
rcorr(as.matrix(data3))

##                          tfr contraception economicActivityFemale
## tfr                     1.00         -0.75                  -0.26
## contraception          -0.75          1.00                   0.15
## economicActivityFemale -0.26          0.15                   1.00
## 
## n= 116 
## 
## 
## P
##                        tfr    contraception economicActivityFemale
## tfr                           0.0000        0.0054                
## contraception          0.0000               0.1190                
## economicActivityFemale 0.0054 0.1190