Regresión y verificación de supuestos

Contexto del conjunto de datos

library(readr)
Fish <- read_csv("Fish.csv")

Estos datos provienen de una muestra de 159 especímenes de 7 especies de peces que son distribuidos en determinado mercado. En estos especímenes se midieron las siguientes variables:

Weight: Peso del pez en gramos
Length 1: Longitud vertical en centímetros
Length 2: Longitud diagonal en centímetros
Length 3: Longitud transversal en centímetros
Height: Altura en centimetros.
Width: Ancho de la diagonal en centímetros.
Species: Especie del pez

El objetivo de este conjunto de datos es explicar el peso del pez a través de las otras variables del conjunto de datos.

El siguiente comando establece un ajuste lineal de los datos usando la variable Weight como dependiente y las demás como explicativas. Se presentan los coeficientes de la regresión.

fit_peces<-lm(Weight~Length1+Length2+Length3+Height+Width+ Species, data=Fish)
coef(fit_peces)

##      (Intercept)          Length1          Length2          Length3 
##      -918.332140       -80.302952        79.888631        32.535381 
##           Height            Width    SpeciesParkki     SpeciesPerch 
##         5.250988        -0.515438       164.722661       137.948910 
##      SpeciesPike     SpeciesRoach     SpeciesSmelt SpeciesWhitefish 
##      -208.429357       103.039955       446.073317        93.874168

De los coeficientes de Species tendremos que estos son puntos de comparación con respecto a los datos de la especie Bream.

anova(fit_peces)

## Analysis of Variance Table
## 
## Response: Weight
##            Df   Sum Sq  Mean Sq   F value    Pr(>F)    
## Length1     1 16978060 16978060 1928.5523 < 2.2e-16 ***
## Length2     1   235134   235134   26.7091 7.589e-07 ***
## Length3     1    74109    74109    8.4181  0.004287 ** 
## Height      1   619028   619028   70.3160 3.787e-14 ***
## Width       1    18474    18474    2.0985  0.149571    
## Species     6  1028534   171422   19.4720 < 2.2e-16 ***
## Residuals 147  1294118     8804                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_peces)$adj.r.squared ## R^2_C

## [1] 0.9313021

Verificamos gráficamente el comportamiento de los residuales

plot(fit_peces)

Verificación de supuestos

Haremos verificación de supuestos básicos en el ajuste del modelo, sobre todo aquellos relacionados con la distribución de la que provienen los datos:

\[e_{ij} \sim^{iid} N(0,\sigma^2)\] Guardamos en un conjunto los valores de los residuos de la regresión

residuos<-resid(fit_peces)

Autocorrelación

En este caso no tiene sentido la aplicabilidad del supuesto de autocorrelación por la forma en que se observaron los especímenes, en caso de contar con un orden intrínseco que tenga relevancia en la obtención de los datos, allí debemos verificar el supuesto, adicionalmente cuando observamos datos con una estructura temporal definida.

El siguiente es el coeficiente de correlación de los residuos.

cor(residuos[-1],residuos[-159])

## [1] 0.5095895

Normalidad

Para verificar el comportamiento distribucional de estos datos debemos tener en cuenta que la estos supuestos deben revisarse a partir del comportamiento de los residuales, los cuales podemos extraerlos directamente desde el modelo.

Se sugiere realizar un gráfico de quantil-quantil para verificar si hay indicios de no normalidad en la muestra seleccionada:

library(car)
qqPlot(residuos, pch = 20)

## [1]  73 143

Hay valores que se salen de las bandas de confianza, por lo tanto, hacemos el cálculo de una prueba de normalidad. Si p-value es menor a 0.05 (dependiendo del caso) decimos que hay evidencia para afirmar que los datos no provienen de una distribución normal.

shapiro.test(residuos)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.93335, p-value = 9.145e-07

Homocedasticidad

Verificamos que la varianza de los datos es igual en todas las particiones definidas para el tipo de especie.

boxplot(residuos~Fish$Species)

bartlett.test(residuos~Fish$Species)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuos by Fish$Species
## Bartlett's K-squared = 16.183, df = 6, p-value = 0.01281

Multicolinealidad

Verificamos sobre las variables independientes si hay alguna que se correlacione directamente con otra, de tal manera que ambas otorguen la misma información.

plot(Fish) # Gráfico de dispersión multivariado

cor.test(Fish$Length1,Fish$Length2)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length1 and Fish$Length2
## t = 403.11, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9993394 0.9996473
## sample estimates:
##       cor 
## 0.9995173

cor.test(Fish$Length1,Fish$Length3)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length1 and Fish$Length3
## t = 98.656, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9891090 0.9941713
## sample estimates:
##      cor 
## 0.992031

cor.test(Fish$Length1,Fish$Height)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length1 and Fish$Height
## t = 10.042, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5203840 0.7117451
## sample estimates:
##       cor 
## 0.6253779

cor.test(Fish$Length1,Fish$Width)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length1 and Fish$Width
## t = 21.806, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8223868 0.9010912
## sample estimates:
##       cor 
## 0.8670497

cor.test(Fish$Length2,Fish$Length3)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length2 and Fish$Length3
## t = 114.86, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9919372 0.9956878
## sample estimates:
##       cor 
## 0.9941026

cor.test(Fish$Length2,Fish$Height)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length2 and Fish$Height
## t = 10.449, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5384698 0.7239251
## sample estimates:
##       cor 
## 0.6404408

cor.test(Fish$Length2,Fish$Width)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length2 and Fish$Width
## t = 22.487, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8308686 0.9060084
## sample estimates:
##       cor 
## 0.8735467

cor.test(Fish$Length3,Fish$Height)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length3 and Fish$Height
## t = 12.4, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6151058 0.7742847
## sample estimates:
##       cor 
## 0.7034089

cor.test(Fish$Length3,Fish$Width)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Length3 and Fish$Width
## t = 23.043, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8373749 0.9097666
## sample estimates:
##       cor 
## 0.8785202

cor.test(Fish$Height,Fish$Width)

## 
##  Pearson's product-moment correlation
## 
## data:  Fish$Height and Fish$Width
## t = 16.303, df = 157, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7269460 0.8443297
## sample estimates:
##      cor 
## 0.792881

Veamos un modelo solo con la variable Length1 y Species.

fit_2<-lm(Weight~Length1+Species,data=Fish)
anova(fit_2)

## Analysis of Variance Table
## 
## Response: Weight
##            Df   Sum Sq  Mean Sq  F value    Pr(>F)    
## Length1     1 16978060 16978060 1810.731 < 2.2e-16 ***
## Species     6  1853569   308928   32.948 < 2.2e-16 ***
## Residuals 151  1415830     9376                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_2)$adj.r.squared

## [1] 0.9268321

plot(fit_2)