Análisis Estadístico con R - Parte II

Requisitos

1. La muestra de datos pareados \((x, y)\) es una muestra aleatoria simple de datos cuantitativos.

2. El examen visual del diagrama de dispersión debe confirmar que los puntos se acercan al patrón de una línea recta.

3. Como los resultados se pueden ver muy afectados por la presencia de valores atípicos, es necesario eliminar cualquier valor atípico, si se sabe que se trata de un error. Los efectos de cualquier otro valor atípico deben tomarse en cuenta calculando \(r\) con y sin el valor atípico incluido.

4. Los pares de datos \((x, y)\) tienen una distribución normal bivariada, este supuesto requiere que, para cualquier valor fijo de \(x\), los valores correspondientes de \(y\) tengan una distribución aproximadamente normal, y que para cualquier valor fijo de \(y\), los valores de \(x\) tengan también una distribución aproximadamente normal.

args(cor)

## function (x, y = NULL, use = "everything", method = c("pearson", 
##     "kendall", "spearman")) 
## NULL

Propiedades

1. El valor de \(r\) está siempre entre \(-1\) y \(1\), inclusive. Es decir \[-1 \le r \le 1\]

2. El valor de \(r\) no cambia si todos los valores de cualquiera de las variables se convierten a una escala diferente.

3. El valor de \(r\) no se ve afectado por la elección de \(x\) o \(y\). Intercambie todos los valores de \(x\) y \(y\), y el valor de \(r\) no sufrirá cambios.

4. \(r\) mide la fuerza de una relación lineal. No está diseñada para medir la fuerza de una relación que no sea lineal.

5. \(r\) es muy sensible a los valores atípicos, en el sentido de que un solo valor atípico puede afectar su valor de manera drástica.

Ejemplo 1

Suponga el siguiente conjunto de datos.

datasets::anscombe

Para los cuatro conjuntos de datos se cumple:

¿Los cuatro conjuntos de datos se comportan igual?

ff <- y ~ x
 for(i in 1:4) {
   ff[2:3] <- lapply(paste(c("y","x"), i, sep=""), as.name)
   assign(paste("lm.",i,sep=""), lmi <- lm(ff, data= anscombe))
 }
 
 op <- par(mfrow=c(2,2), mar=1.5+c(4,3.5,0,1), oma=c(0,0,0,0),
           lab=c(6,6,7), cex.lab=1.5, cex.axis=1.3, mgp=c(3,1,0))

 for(i in 1:4) {
   ff[2:3] <- lapply(paste(c("y","x"), i, sep=""), as.name)
   plot(ff, data =anscombe, col="red", pch=21, bg = "orange", cex = 2.5,
        xlim=c(3,19), ylim=c(3,13),las=1,
        xlab=eval(substitute(expression(x[i]), list(i=i))),
        ylab=eval(substitute(expression(y[i]), list(i=i))))
   abline(get(paste("lm.",i,sep="")), col="blue")
 }

Ejemplo 2

if(!require(datasauRus)) install.packages("datasauRus",dependencies = T)
library(dplyr)
datasaurus_dozen%>% 
    group_by(dataset) %>% 
    summarize(
      mean_x    = mean(x),
      mean_y    = mean(y),
      std_dev_x = sd(x),
      std_dev_y = sd(y),
      corr_x_y  = cor(x, y)
    )

require(plotly)
datasaurus_dozen %>%
  plot_ly(x=~x, y=~y,
          color=~dataset,
          frame=~dataset,
          type="scatter",
          mode="markers") %>%
  animation_opts(2000,easing = "exp", redraw = FALSE)

Errores comunes en la interpretación de la correlación

1. Un error común es concluir que la correlación implica causalidad. Si el aumento de los conciertos de Madonna en los años 80 hace que disminuya la población de ballenas en el pacífico, ¿Significa que Madonna es una asesina de ballenas? En algunos casos, los datos pueden verse afectados por alguna otra variable interventora en los antecedentes. (Una variable interventora es aquella que afecta a las variables que se estudian, pero que no está incluida en la investigación).

2. Otro error proviene de los datos basados en promedios. Los promedios eliminan la variación individual y pueden inflar el coeficiente de correlación. Un estudio produjo un coeficiente de correlación lineal de 0.4 para datos pareados que relacionaban el ingreso y la escolaridad de individuos, pero el coeficiente de correlación lineal se convirtió en 0.7 cuando se utilizaron promedios regionales.

3. Un tercer error implica la propiedad de linealidad. Si no hay una correlación lineal, podría existir una correlación no lineal.

Requisitos

1. La muestra de datos pareados \((x, y)\) es una muestra aleatoria de datos cuantitativos.

2. El examen visual del diagrama de dispersión indica que los puntos se aproximan al patrón de una línea recta.

3. Los valores atípicos pueden tener un gran efecto sobre la ecuación de regresión, por lo que se debe eliminar cualquier valor atípico, si se sabe que es un error. Es importante tomar en cuenta los efectos de cualquier valor atípico que no sea un error conocido.

4. Para cada valor fijo de \(x\), los valores correspondientes de \(y\) tienen una distribución en forma de campana.

5. Para los distintos valores fijos de \(x\), las distribuciones de los valores correspondientes de \(y\) tienen la misma desviación estándar.

6. Para los distintos valores fijos de \(x\), las distribuciones de los valores correspondientes de \(y\) tienen medias que se ubican en la misma línea recta.

Supuestos

• Linealidad: la relación entre \(X\) e \(Y\) debe ser lineal.

• Homoscedasticidad: la varianza de los residuales \(\varepsilon\) son iguales para cualquier valor de \(X\).

• Independencia: todas las observaciones son independientes unas de otras.

• Normalidad: La variable \(Y\) está normalmente distribuída para cualquier valor de \(X\).

\[y \sim \beta_0 + \beta_1x+\varepsilon \quad\text{donde}\quad \varepsilon\sim N(0,\sigma_\varepsilon^2)\]

Residuales

El criterio utilizado para determinar cuál recta es “mejor”" que todas las demás se basa en las distancias verticales entre los puntos de datos originales y la recta de regresión. Tales distancias se denominan residuos.

Para una muestra de datos pareados \(x\) y \(y\), un residuo es la diferencia entre un valor y muestral observado y el valor de y predicho por medio de la ecuación de regresión. Es decir, \[residuo=y_\text{observada}- y_\text{predicha}=y-\hat{y}\]

Ejemplo

require(ggplot2)
mtcars

mtcars %>%
  ggplot(aes(x=wt,y=mpg))+geom_point(size=2)+
  xlab("Peso del Vehículo (1000 lbs)")+
  ylab("Millas/Galón(US)")+
  ggtitle("Relación entre el peso de los vehículos y el recorrido")+
  theme_bw()

with(mtcars,cor(y = mpg,x = wt))

## [1] -0.8676594

with(mtcars,cor.test(y = mpg,x = wt))

## 
##  Pearson's product-moment correlation
## 
## data:  wt and mpg
## t = -9.559, df = 30, p-value = 1.294e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9338264 -0.7440872
## sample estimates:
##        cor 
## -0.8676594

args(lm)

## function (formula, data, subset, weights, na.action, method = "qr", 
##     model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, 
##     contrasts = NULL, offset, ...) 
## NULL

fit <- lm(mpg~wt,data = mtcars)
fit

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Coefficients:
## (Intercept)           wt  
##      37.285       -5.344

coef(fit)

## (Intercept)          wt 
##   37.285126   -5.344472

mtcars %>%
  ggplot(aes(x=wt,y=mpg))+
  geom_abline(intercept = coef(fit)[1],
              slope = coef(fit)[2],
              colour=2,size=1.5)+
  geom_point(size=2)+
  xlab("Peso del Vehículo (1000 lbs)")+
  ylab("Millas/Galón(US)")+
  ggtitle("Relación entre el peso de los vehículos y el recorrido")+
  theme_bw()

(res <- residuals(fit))

##           Mazda RX4       Mazda RX4 Wag          Datsun 710      Hornet 4 Drive 
##          -2.2826106          -0.9197704          -2.0859521           1.2973499 
##   Hornet Sportabout             Valiant          Duster 360           Merc 240D 
##          -0.2001440          -0.6932545          -3.9053627           4.1637381 
##            Merc 230            Merc 280           Merc 280C          Merc 450SE 
##           2.3499593           0.2998560          -1.1001440           0.8668731 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood Lincoln Continental 
##          -0.0502472          -1.8830236           1.1733496           2.1032876 
##   Chrysler Imperial            Fiat 128         Honda Civic      Toyota Corolla 
##           5.9810744           6.8727113           1.7461954           6.4219792 
##       Toyota Corona    Dodge Challenger         AMC Javelin          Camaro Z28 
##          -2.6110037          -2.9725862          -3.7268663          -3.4623553 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2        Lotus Europa 
##           2.4643670           0.3564263           0.1520430           1.2010593 
##      Ford Pantera L        Ferrari Dino       Maserati Bora          Volvo 142E 
##          -4.5431513          -2.7809399          -3.2053627          -1.0274952

cbind(mean(res),var(res))

##               [,1]     [,2]
## [1,] -5.114724e-17 8.978127

plot(res,pch=19,xlab = "Observación",
     ylab="Residuales del Modelo",las=1,
     main="Comportamiento de los residuales del modelo")
abline(h=0,col=2)

hist(res,col="seagreen3",freq = F,las=1,
     xlab="Residuales del modelo",ylab="Densidad",
     main="Comportamiento de los residuales",
     xlim=c(-10,10))
curve(dnorm(x,mean = mean(res),sd = sd(res)),col=4,lwd=2,add=T)

plot(fit,pch=19,cex=1.2)

shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.94508, p-value = 0.1044

if(!require(car)) install.packages("car",dependencies = T)
outlierTest(fit)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##          rstudent unadjusted p-value Bonferroni p
## Fiat 128 2.537801           0.016788       0.5372

leveragePlots(fit,pch=19)

# Gráfico  Distacia Cook's, identifique valores D > 4/(n-k-1) 
(cutoff <- 4/((nrow(mtcars)-length(fit$coefficients)-1)) )

## [1] 0.137931

plot(fit, which=4, cook.levels=cutoff)

influencePlot(fit,id.method="identify",
              main="Gráfico de Influencia",
              sub="El tamaño del círcuo es la proporción de la distancia Cook's")

ncvTest(fit)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.03794177, Df = 1, p = 0.84556

spreadLevelPlot(fit,pch=19)

## 
## Suggested power transformation:  1.347793

crPlots(fit,pch=19)

durbinWatsonTest(fit)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.3628798      1.251727   0.028
##  Alternative hypothesis: rho != 0

if(!require(gvlma)) install.packages("gvlma",dependencies = T)
gvmodel <- gvlma(fit) 
summary(gvmodel)

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = fit) 
## 
##                       Value  p-value                   Decision
## Global Stat        11.73816 0.019408 Assumptions NOT satisfied!
## Skewness            2.37864 0.123004    Assumptions acceptable.
## Kurtosis            0.02033 0.886622    Assumptions acceptable.
## Link Function       8.57441 0.003409 Assumptions NOT satisfied!
## Heteroscedasticity  0.76478 0.381838    Assumptions acceptable.

Estimación de los parámetros

Al igual que en el caso de la regresión lineal simple, el objetivo de la estimación de mínimos cuadrados es encontrar valores para los parámetros \(\beta_0,\beta_1,...,\beta_p\) que minimicen el error estándar. Este valor se denotará como \(Q\) y se define como:

\[Q=\sum_{i=1}^{n}(Y_i-\beta_0-\beta_1x_1-\dots-\beta_px_p)^2\]

El vector (de los parámetros estimados se denota como:

\[\underset{p\times 1}{b}=\begin{bmatrix} b_0\\ b_1\\ \vdots\\ b_p\\ \end{bmatrix}\]

Por lo tanto, las ecuaciones de mínimos cuadrados para el modelo de regresión múltiple se pueden representar como:

\[X^\prime Xb=X^\prime Y\]

Con estimadores

\[\underset{2\times 1}{b}=(\underset{2\times 2}{X^\prime X})^{-1}(\underset{2\times 1}{X^\prime X})Y\]

La respuesta ajustada también se puede representar como un vector denotado \(\hat{Y}\) con valores \(\hat{Y_i}\) y residuos \(\varepsilon_i=Y_i − \hat{Y_i}\).

\[\underset{n\times 1}{\hat{Y}}=\begin{bmatrix} \hat{Y_1}\\ \hat{Y_2}\\ \vdots\\ \hat{Y_n}\\ \end{bmatrix} \quad \underset{n\times 1}{\varepsilon}=\begin{bmatrix} \varepsilon_1\\ \varepsilon_2\\ \vdots\\ \varepsilon_n\\ \end{bmatrix}\]

Los valores ajustados también se pueden escribir como:

\[\underset{n\times 1}{\hat{Y}}=Xb\]

Con residuales

\[\underset{n\times 1}{\varepsilon}=Y-\hat{Y}=Y-Yb\]

La matriz gorro es definida como

\[H=X(X^\prime X)^{-1}X^\prime\]

Por lo tanto, el vector de valores ajustados \(\hat{Y}\) se puede expresar en términos de la matriz gorro:

\[\underset{n\times 1}{\hat{Y}}=HY\]

Con vector de residuales

\[\underset{n\times 1}{\varepsilon}=(I-H)Y\]

Análisis de Varianza (ANOVA)

Al referirse brevemente ANOVA en la configuración de regresión múltiple, la suma de los cuadrados y la media de los cuadrados de cada fuente de variación se definen en términos de matriz como:

\[\begin{align*} SSE = & Y^\prime Y-b^\prime X^\prime Y=Y^\prime(I-H)Y\\ SST = & Y^\prime Y-\frac{\sum(Y^2)}{n}\\ SSR = & SST-SSE\\ MSE = & \frac{SSE}{n-p}\\ MSR = & \frac{SSR}{p-1} \end{align*}\]

Ejemplo

library(car)
library(MASS)
library(gvlma)
library(ggplot2)
library(gridExtra)
mtcars$am<-as.factor(mtcars$am)
levels(mtcars$am)<-c('Automatic','Manual')
tapply(mtcars$mpg,mtcars$am,mean)

## Automatic    Manual 
##  17.14737  24.39231

mtcars %>% 
  ggplot(aes(x=am,y=mpg,fill=am))+geom_boxplot()

var.test(mpg~am,data = mtcars)

## 
##  F test to compare two variances
## 
## data:  mpg by am
## F = 0.38656, num df = 18, denom df = 12, p-value = 0.06691
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1243721 1.0703429
## sample estimates:
## ratio of variances 
##          0.3865615

t.test(mpg~am,data=mtcars,var.equal=T)

## 
##  Two Sample t-test
## 
## data:  mpg by am
## t = -4.1061, df = 30, p-value = 0.000285
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -10.84837  -3.64151
## sample estimates:
## mean in group Automatic    mean in group Manual 
##                17.14737                24.39231

fit<-lm(mpg~.,data=mtcars)
summary(fit)$coef

##                Estimate  Std. Error    t value   Pr(>|t|)
## (Intercept) 12.30337416 18.71788443  0.6573058 0.51812440
## cyl         -0.11144048  1.04502336 -0.1066392 0.91608738
## disp         0.01333524  0.01785750  0.7467585 0.46348865
## hp          -0.02148212  0.02176858 -0.9868407 0.33495531
## drat         0.78711097  1.63537307  0.4813036 0.63527790
## wt          -3.71530393  1.89441430 -1.9611887 0.06325215
## qsec         0.82104075  0.73084480  1.1234133 0.27394127
## vs           0.31776281  2.10450861  0.1509915 0.88142347
## amManual     2.52022689  2.05665055  1.2254035 0.23398971
## gear         0.65541302  1.49325996  0.4389142 0.66520643
## carb        -0.19941925  0.82875250 -0.2406258 0.81217871

round(vif(fit),2)

##   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb 
## 15.37 21.62  9.83  3.37 15.16  7.53  4.97  4.65  5.36  7.91

sqrt(vif(fit)) > 2

##   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb 
##  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

stepAIC(fit,direction='both',trace=F,scope=list(upper=~.,lower=~1))$anova

BestFit<-update(fit,.~wt+qsec+am)
summary(BestFit)$coef

##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  9.617781  6.9595930  1.381946 1.779152e-01
## wt          -3.916504  0.7112016 -5.506882 6.952711e-06
## qsec         1.225886  0.2886696  4.246676 2.161737e-04
## amManual     2.935837  1.4109045  2.080819 4.671551e-02

vif(BestFit)

##       wt     qsec       am 
## 2.482952 1.364339 2.541437

sqrt(vif(BestFit)) > 2

##    wt  qsec    am 
## FALSE FALSE FALSE

FinalFit<-update(fit,.~am:wt+am:qsec)
vif(FinalFit)

##             GVIF Df GVIF^(1/(2*Df))
## am:wt   17.50104  2        2.045342
## am:qsec 17.50104  2        2.045342

spreadLevelPlot(FinalFit,pch=19)

## 
## Suggested power transformation:  1.035865

par(mfrow=c(1,2))
qqnorm(resid(FinalFit),pch=19);qqline(resid(FinalFit))
plot(fitted(FinalFit),resid(FinalFit),pch=19,
     main='Residuals vs Fitted',xlab='Fitted',ylab='Residuals')
abline(h=0,col='red')

leveragePlots(FinalFit,pch=19)

avPlots(FinalFit,pch=19)

ncvTest(FinalFit)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.6460701, Df = 1, p = 0.42152

durbinWatsonTest(FinalFit)

##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1114435      2.164208   0.952
##  Alternative hypothesis: rho != 0

summary(gvlma(FinalFit))

## 
## Call:
## lm(formula = mpg ~ am:wt + am:qsec, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9361 -1.4017 -0.1551  1.2695  3.8862 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       13.9692     5.7756   2.419  0.02259 *  
## amAutomatic:wt    -3.1759     0.6362  -4.992 3.11e-05 ***
## amManual:wt       -6.0992     0.9685  -6.297 9.70e-07 ***
## amAutomatic:qsec   0.8338     0.2602   3.205  0.00346 ** 
## amManual:qsec      1.4464     0.2692   5.373 1.12e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.097 on 27 degrees of freedom
## Multiple R-squared:  0.8946, Adjusted R-squared:  0.879 
## F-statistic: 57.28 on 4 and 27 DF,  p-value: 8.424e-13
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = FinalFit) 
## 
##                      Value p-value                Decision
## Global Stat        3.13417  0.5356 Assumptions acceptable.
## Skewness           0.51648  0.4723 Assumptions acceptable.
## Kurtosis           0.48168  0.4877 Assumptions acceptable.
## Link Function      2.07096  0.1501 Assumptions acceptable.
## Heteroscedasticity 0.06506  0.7987 Assumptions acceptable.

Anova(FinalFit)

p1 <- ggplot(mtcars,aes(x=wt,y=mpg,colour=am))+
  geom_point()+geom_smooth(method = "lm",se = F)+
  theme(legend.position="bottom")
p2 <- ggplot(mtcars,aes(x=qsec,y=mpg,colour=am))+
  geom_point()+geom_smooth(method = "lm",se = F)

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mylegend<-g_legend(p1)

grid.arrange(arrangeGrob(p1 + theme(legend.position="none"),
                         p2 + theme(legend.position="none"),
                         nrow=1),
             mylegend, nrow=2,heights=c(10, 1))