Lugo_Ferran_Reg_PEC1
Ferran Lugo
2023-04-22
PEC 1: REGRESIÓN, MODELOS Y MÉTODOS
Alumno: Ferran Lugo UOC - Máster en Bioinformática y Bioestadística
Ejercicio 1
Un grupo de científicos norteamericanos están interesados en encontrar un hábitat adecuado para reintroducir una especie rara de escarabajos tigre, llamada cicindela dorsalis dorsalis, los cuales viven en playas de arena de la costa del Atlántico Norte. Se muestrearon 12 playas y se midió la densidad de estos escarabajos tigre. Adicionalmente se midieron una serie de factores bióticos y abióticos tales como la exposición a las olas, tamaño de la partícula de arena, pendiente de la playa y densidad de los anfípodos depredadores. Los datos se hallan en la hoja de cálculo cicindela.xlsx.
(a) Ajustar un modelo de regresión lineal múltiple que estime todos los coeficientes de regresión parciales referentes a todas las variables regresoras y el intercepto.
Primeramente, cargaremos los datos y definiremos las variables.
library(readxl)## Warning: package 'readxl' was built under R version 4.2.3cicindela <- read_excel("C:/Users/ferra/OneDrive/Escritorio/cicindela.xlsx")Las variables de trabajo son:
BeetleDensity: Densidad de escarabajos (no se especifica unidades).
Wave exposure: Exposición a las olas (no se especifica unidades)
Sandparticlesize: Tamaño de partícula de arena (no se especifica unidades, probablemente mm)
Beachsteepness: Pendiente de la playa (no se especifica unidades, probablemente grados)
AmphipodDensity: Densidad de depredadores (no se especifica unidades)
Según el enunciado, la variable que se desea analizar respecto las demás (variable dependiente) será BeetleDensity. y todas las demás se incluirán a priori como variables independientes. Así, como existen cuatro variables independientes, el modelo de regresión lineal sería:
Y = beta1 * X1 + beta2 * X2 + beta3 * X3 + beta4 * X4 + ei
Donde Y representa la variable dependiente (BeetleDensity); las “betas” son el efecto promedio que tiene el incremento en una unidad de la variable independiente en cuestión sobre la dependiente, las “X” son las respectivas variables independientes (por orden, Wave exposure, Sandparticlesize, BeachSteepness y AmphipodDensity) y “e” representa el posible error existente.
¿Es significativo el modelo obtenido? ¿Qué test estadístico se emplea para contestar a esta pregunta. Plantear la hipótesis nula y la alternativa del test.
Así, crearemos el modelo de regresión con comandos (lm).
# Mediante "." hacemos incluir el resto de variables del dataframe como independientes.
beetle_lm <- lm(BeetleDensity ~ ., data = cicindela)
summary(beetle_lm)##
## Call:
## lm(formula = BeetleDensity ~ ., data = cicindela)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3004 -2.7038 0.0795 2.6017 5.3924
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.9531 17.2661 0.866 0.4152
## `Wave exposure` 0.9123 1.0935 0.834 0.4317
## Sandparticlesize 3.8970 1.1690 3.334 0.0125 *
## `Beach steepness` 0.6511 0.4530 1.437 0.1938
## AmphipodDensity -1.5624 0.6610 -2.364 0.0501 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.513 on 7 degrees of freedom
## Multiple R-squared: 0.9578, Adjusted R-squared: 0.9337
## F-statistic: 39.71 on 4 and 7 DF, p-value: 6.727e-05La significación del modelo la podemos analizar mediante el p-valor asociado al test F del modelo (última línea), que indica la significación estadística de las relaciones que describe el modelo. Este pvalor permite determinar cuál de las dos hipótesis elegir en base al siguiente contraste:
Ho (hipótesis nula): Las betas del modelo son iguales cero
H1 (hipótesis alternativa): Como mínimo alguna de las betas no es igual cero
Recordemos que este pvalor se obtiene mediante el test F, cuyo estadístico en este caso es 39.71. Éste consiste en una relación de la varianza explicada por nuestro modelo de regresión en relación con un modelo con solo el intercepto y ninguna otra variable. En este caso, como el p-valor es inferior a 0.05 (que es el nivel de significación por defecto, alfa), rechazaremos la hipótesis nula (Ho), aceptando la hipótesis alternativa (H1). Es decir, alguna de las betas del modelo no es cero y, por lo tanto, las relaciones establecidas en el modelo pueden ser significativas.
En la primera columna observamos todos los coeficientes de regresión para cada variable, más el del intercepto (el primero). En este caso, la única relación significativa se encontraría entre la variable dependiente (BeetleDensity) y Sandparticlesize, siendo una relación positiva con un estimador de 3.8970. En otras palabras, por cada aumento de una unidad de la medida de la arena se espera un aumento de 3.897 unidades de densidad de escarabajos.
Hace falta remarcar que se está siendo estricto, puesto que existe un p-valor de 0.0501 entre AmphipodDensity y BeetleDensity que posiblemente pueda ser importante si se cambia el nivel de significación alfa.
Otra variable interesante a analizar sería el “Multiple R-squared”, que identifica la proporción de variación explicada por este modelo. Al ser muy próxima a 1 (0.9578), podríamos aceptar que este modelo se ajusta a los datos. En otras palabras, un 95.78% de la variación de la variable respuesta puede ser explicada por las variables predictoras.
¿Qué variables han salido significativas para un nivel de significación α = 0.10?
En el análisis anterior se nos muestan los p-valores que marcan la significación de la relación para cada variable (columna de la derecha). Si tomamos un nivel de significación de alfa = 0.1, deberemos analizar aquellas variables con p-valores inferiores a 0.1 en el “summary”. Así, las variables significativas para este nivel de significación serían Sandparticlesize (p-valor 0.0125) y AmphipodDensity (p-valor 0.0501). El resto tiene p-valores de más de 0.1, por lo que no entrarían dentro de nuestro baremo.
(b) Calcular los intervalos de confianza al 90 y 95% para el parámetro que acompaña a la variable AmphipodDensity. Utilizando sólo estos intervalos, ¿qué podríamos haber deducido sobre el pvalor para la densidad de los anfípodos depredadores en el resumen del modelo de regresión? ¿Qué interpretación práctica tiene este parámetro β4?
Mediante el comando confint podemos calcular los intervalos de confianza para todos los parámetros de las variables del modelo.
# Al 90%
confint(beetle_lm,level=.90)## 5 % 95 %
## (Intercept) -17.7588417 47.6650535
## `Wave exposure` -1.1594063 2.9840266
## Sandparticlesize 1.6823301 6.1117347
## `Beach steepness` -0.2070857 1.5093046
## AmphipodDensity -2.8146991 -0.3100058# Al 95%
confint(beetle_lm,level=.95)## 2.5 % 97.5 %
## (Intercept) -25.8746879 5.578090e+01
## `Wave exposure` -1.6733999 3.498020e+00
## Sandparticlesize 1.1328616 6.661203e+00
## `Beach steepness` -0.4200042 1.722223e+00
## AmphipodDensity -3.1254068 7.019125e-04Así, los intervalos al 90% para la variable AmphipodDensity son: (-2.8146991, -0.3100058).
Y los intervalos para la misma variable al 95%, serían: (-3.1254068, 7.019125e-04).
Debido a que el intervalo de confianza al 95% incluye el valor cero, pero el intervalo al 90%, no, podemos deducir que el parámetro de la variable en cuestión (AmphipodDensity) es signigicativamente distinta a cero para la significación del 10% (primer caso), pero no para la significación del 5% (segundo caso). Este hecho es coherente con el resultado del “summary” anterior, donde hemos visto que para el nivel de significación alfa = 0.1 esta variable era significativa, pero no para el nivel de significación alfa = 0.05.
Este resultado nos permite verificar que la densidad de depredadores (AmphipodDensity) puede tener un efecto sobre la variable regresora. Si nos fijamos en su valor estimado en el “summary” (beta4 = -1.5624), significaría que, para cada unidad de adición de densidad de depredadores (AmphipodDensity), se espera un decremento de -1.5624 unidades de la variable dependiente, la densidad de escarabajos (BeetleDensity). Este resultado tiene lógica, debido a que los Amfípodos son depredadores de los escarabajos (cuanto más depredadores, menos presas habrá).
(c) Estudiar la posible multicolinealidad del modelo con todas las regresoras calculando los VIFs.
Como tomaremos en cuenta todas las variables, seguimos con el modelo anterior. Calcularemos los VIF para cada variable regresora del modelo, gracias al paquete car y a la función vif. Los representaremos en un barplot, especificando el umbral de multicolinealidad en valores de vif de 5 o más (es lo recomendado en la bibliografía).
library(car)## Warning: package 'car' was built under R version 4.2.2## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.2vifs_beetle <- vif(beetle_lm)
vifs_beetle## `Wave exposure` Sandparticlesize `Beach steepness` AmphipodDensity
## 3.771652 3.398998 1.158425 5.119632# Los representamos en un bar plot
barplot(vifs_beetle, col = "lightblue", horiz = TRUE)
abline(v = 5, lwd = 3, lty = 2)
El VIF para Wave exposure es 3.771652; para Sandparticlesize, 3.398998; para Beach steepness, 1.158425 y, finalmente, para AmphipodDensity, 5.119632.
La multicolinealidad ocurre cuando dos o más variables están altamente correlacionadas, afectando a la viabilidad del modelo (ya que por defecto se relacionan con la variable independiente por separado). Calculando los VIFs podemos averiguar si existe tal caso. Un VIF de 1 indica que no existe correlación alguna entre la variable predictora en cuestión y otra variable predictora del modelo, mientras que un VIF mayor que 5 indica que existe potencialmente una correlación entre la variable en cuestión y otra variable del modelo.
En este caso, observamos que AmphipodDensity tiene un VIF superior a 5, sugeriendo que es probable que esté correlacionada con otra variable del modelo. Para detectar con qué variable puede estar correlacionada, podemos analizar la matriz de correlaciones de las variables en el dataframe inicial.
cor(cicindela)## BeetleDensity Wave exposure Sandparticlesize Beach steepness
## BeetleDensity 1.0000000 0.84898795 0.91454493 0.21821349
## Wave exposure 0.8489880 1.00000000 0.77117855 0.05235125
## Sandparticlesize 0.9145449 0.77117855 1.00000000 0.01310056
## Beach steepness 0.2182135 0.05235125 0.01310056 1.00000000
## AmphipodDensity -0.9348247 -0.83984838 -0.81540614 -0.20522682
## AmphipodDensity
## BeetleDensity -0.9348247
## Wave exposure -0.8398484
## Sandparticlesize -0.8154061
## Beach steepness -0.2052268
## AmphipodDensity 1.0000000# Lo representamos con un mapa de calor
heatmap(cor(cicindela))
Observamos que la variable AmphipodDensity tiene una correlación negativa muy elevada con BeetleDensity (-0.9348247, muy próxima a -1), cuyo marcaje en el mapa de calor indica que pueden ser las variables con multicolinealidad en este modelo. Para solucionar la multicolinealidad, podríamos plantearnos quitar la variable problemática del modelo.
(d) Considerar el modelo más reducido que no incluye las variables exposición a las olas y la pendiente de la playa y decidir si nos podemos quedar con este modelo reducido mediante un contraste de modelos con el test F para un α = 0.05.
Primeramente, creamos el modelo del enunciado y anotamos el summary para observar su significación individual.
beetle_lm2 <- lm(BeetleDensity ~ Sandparticlesize + AmphipodDensity, data = cicindela)
summary(beetle_lm2)##
## Call:
## lm(formula = BeetleDensity ~ Sandparticlesize + AmphipodDensity,
## data = cicindela)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.933 -2.226 -0.512 3.315 5.787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.5651 9.4259 3.773 0.00440 **
## Sandparticlesize 3.7103 1.1215 3.308 0.00911 **
## AmphipodDensity -2.1228 0.5167 -4.108 0.00264 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.621 on 9 degrees of freedom
## Multiple R-squared: 0.9431, Adjusted R-squared: 0.9305
## F-statistic: 74.58 on 2 and 9 DF, p-value: 2.501e-06Observamos que, individualmente, este modelo tiene un p-valor inferior al nivel de significación alfa=0.05 y, por lo tanto, puede ser explicativo.
Pero como nos interesa el contraste entre modelos, los compararemos con el test F. Para ello utilizaremos el comando anova con ambos. La hipótesis nula (Ho) para éste será que los parámetros no incluidos en el primer modelo son cero, mientras que la hipótesis alternativa (H1) será que almenos uno de ellos no lo es.
anova(beetle_lm,beetle_lm2)## Analysis of Variance Table
##
## Model 1: BeetleDensity ~ `Wave exposure` + Sandparticlesize + `Beach steepness` +
## AmphipodDensity
## Model 2: BeetleDensity ~ Sandparticlesize + AmphipodDensity
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 7 142.59
## 2 9 192.19 -2 -49.61 1.2178 0.3517Con un p-valor de 0.3517, superior al nivel de significación alfa = 0.05, aceptamos la hipótesis nula (Ho) y, por lo tanto, podemos afirmar que los parámetros de las variables no incluídas en el modelo simple son cero. Así, parece ser que el modelo simplificado puede ser utilizable.
Escribir en forma paramétrica las hipótesis H0 y H1 de este contraste. Comparar el ajuste de ambos modelos.
Como se ha mencionado, la forma paramétrica incluiría que los parámetros que hemos eliminado en el segundo modelo son cero como hipótesis nula (Ho), mientras que estos parámetros serán distintos a cero para la hipótesis alternativa (H1). Como las variables que hemos eliminado son Wave exposure (primera variable independiente del primer modelo) y Beach steepness (tercera variable independiente del primer modelo), especificaremos sus betas:
Ho: β1 = β3 = 0 H1: β1 o β3 ≠ 0
La comparación del ajuste de ambos modelos la podemos observar en el resumen de cada uno de ellos.
Como hemos visto, el primer modelo (con todas las variables), tenía un p-valor de 6.727e-05 y una R^2 de 0.9578, mientras que este segundo modelo más acotado tiene un p-valor de 2.501e-06, y su R^2 es ligeramente inferior, 0.9431.
Por lo tanto, aunque ambos modelos son significativos (tienen p-valores inferiores al nivel de significación estándar alfa = 0.05), parece que el primer modelo está mejor ajustado que el segundo, al tener un multiple R-squared superior y por lo tanto más próxima a 1. Como este parámetro indica la proporción de la varianza en la variable de respuesta que puede ser explicada por las variables predictoras, parece que el primer modelo sigue siendo el más explicativo.
(e) Calcular y dibujar una región de confianza conjunta al 95% para los parámetros asociados con Sandparticlesize y AmphipodDensity con el modelo que resulta del apartado anterior. Dibujar el origen de coordenadas. La ubicación del origen respecto a la región de confianza nos indica el resultado de una determinada prueba de hipótesis. Enunciar dicha prueba y su resultado.
Calcularemos la región de confianza gracias al paquete ellipse, especificando las variables del modelo. Tal y como están presentadas, se encuentran en las posiciones 2 y 3 del modelo.
require(ellipse)## Loading required package: ellipse## Warning: package 'ellipse' was built under R version 4.2.2##
## Attaching package: 'ellipse'## The following object is masked from 'package:car':
##
## ellipse## The following object is masked from 'package:graphics':
##
## pairs# Dibujamos la región de confianza para el modelo.
# Como observaremos, el origen no está dentro del modelo, por lo que aumentaremos el abarcamiento de los ejes para poderlo dibujar.
plot(ellipse(beetle_lm2,c(2,3)), type="l", xlim=c(-1,7), ylim=c(-4,1))
# Representamos los coeficientes mediante puntos.
points(coef(beetle_lm2)[2], coef(beetle_lm2)[3])
# Representamos el centro y añadimos una etiqueta arriba del mismo (0,0) para reconocerlo
points(0,0)
text(0,0, labels="(0,0)",pos=3)
Observamos que el origen (0,0) no se encuentra dentro de la región de confianza al 95%. El contraste de hipótesis que generaría este resultado sería:
Hipótesis nula (Ho): el parámetro beta de Sandparticlesize es el mismo que el de AmphipodDensity y ambos son nulos
H0: β(Sandparticlesize) = β(AmphipodDensity) = 0
Hipótesis alternativa (H1): Almenos uno de estos parámetros no es nulo.
H1: β(Sandparticlesize) ≠ 0 y/o β(AmphipodDensity) ≠ 0
Como el origen se encuentra fuera de la región de confianza, aceptamos la hipótesis alternativa y, por lo tanto, parece significativo que almenos uno de los coeficientes beta de estas variables no sea cero.
(f) Con el modelo reducido del apartado (d), predecir en forma de intervalo de confianza al 95% la densidad de los escarabajos tigre previsible para una playa cercana a un conocido hotel donde el tamaño de partícula de arena es 5 y la densidad de anfípodos depredadores es 11. Comprobar previamente que los valores observados no suponen una extrapolación.
Según el enunciado tenemos que: Sandparticlesize = 5 AmphipodDensity = 11
Para demostrar que una sandparticle de 5 y una AmphipodDensity de 11 no implican una predicción extrapolable, debemos fijarnos en si son números que queden fuera de los valores sobre los cuales se ha construido el modelo de ajuste. Para ello nos fijaremos en los máximos y mínimos que incluimos en nuestra muestra para ambas variables.
max(cicindela$Sandparticlesize)## [1] 7min(cicindela$Sandparticlesize)## [1] 1max(cicindela$AmphipodDensity)## [1] 19min(cicindela$AmphipodDensity)## [1] 5Observamos que 5 está entre 1 (mínimo) y 7 (máximo) para la variable “Sandparticle”, mientras que 11 está entre 5 (mínimo) y 19 (máximo) para la variable “AmphipodDensity”. Por lo tanto, la predicción con estos valores no sería una extrapolación.
Creamos la playa del enunciado como un dataframe, y realizamos la predicción de su densidad de escarabajos con el modelo anterior. Al anotar “interval=prediction”, se nos mostrará un intervalo de confianza por defecto.
playa_hotel <- data.frame(Sandparticlesize=5,AmphipodDensity=11)
predict(beetle_lm2,playa_hotel,interval="prediction")## fit lwr upr
## 1 30.76569 19.29834 42.23304Así, con el modelo beetle_lm2, que incluye como variables independientes el grueso de la arena y la densidad de anfípotos, se espera que en esa playa haya desde 19.298 unidades hasta 42.233 unidades de densidad de escarabajos.
Ejercicio 2
En el trabajo de Whitman et al. (2004) se estudia, entre otras cosas, la relación entre la edad de los leones y la proporción oscura en la coloración de sus narices. En el archivo lions.csv disponemos de los datos de 105 leones machos y hembras de dos áreas de Tanzania, el parque nacional de Serengueti y el cráter del Ngorongoro, entre 1999 y 2002. Las variables registradas son la edad conocida de cada animal y la proporción oscura de su nariz a partir de fotografías tratadas digitalmente (ver figura adjunta). En la figura 1 se reproduce el gráfico de dispersión de la figura 4 del artículo con el cambio de coloración de la nariz según la edad de machos y hembras en las dos poblaciones separadas. Nota: Los datos se han extraído principalmente del gráfico del artículo de Whitman et al. (2004) y por lo tanto son aproximados. Algunos paquetes de R contienen un data.frame con una parte de estos datos. Por ejemplo LionNoses del paquete abd contiene los datos de todos los machos. En consecuencia, los resultados numéricos de vuestro análisis pueden ser ligeramente distintos a los del trabajo original.
(a) Reproducir el gráfico de dispersión de la figura 1 (figura 4d del artículo) lo más fielmente posible al original, ya que se trata de una exigencia de los editores de la revista.
Primeramente cargamos el archivo y definimos las variables.
library(readr)## Warning: package 'readr' was built under R version 4.2.1lions <- read_csv("C:/Users/ferra/OneDrive/Escritorio/lions.csv")## Rows: 105 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): sex, area
## dbl (2): prop.black, age
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.prop.black: Denota la proporción total de área negra de la nariz (área en negro/área total).
age: Anota la edad del sujeto (león) en años, probablemente.
sex: Sexo del sujeto (león). M es macho y F, hembra.
area: Área de estudio. S es Serengeti y N, Ngorongoro.
Para la realización del plot, tendremos en cuenta que la variable en el eje X es la edad, mientras que la variable en el eje Y es la proporción de color negro. Para que reconozca dos variables como una sola para la leyenda, juntaremos el area con el sexo de los sujetos en un nuevo dataframe (para no comprometer el dataframe inicial).
# Creamos un nuevo dataframe que es igual que el inicial
lions_merge <- lions
# Juntamos las dos últimas columnas
lions_merge$group <- paste(lions_merge$sex,lions_merge$area)
# Quitamos las variables redundantes
lions_merge <- subset(lions_merge, select=c("prop.black","group","age"))En base a estos nuevos cuatro grupos (fusiones de sexo-área), delimitaremos las formas de cada grupo con “pch”, según la bibliografía de formas encontrada en “http://www.sthda.com/english/wiki/r-plot-pch-symbols-the-different-point-shapes-available-in-r” (serán las referencias 1, 2, 16 y 17, es decir, redondas y triángulos blancos/negros). Seguidamente, añadiremos la leyenda abajo a la derecha, anotando los títulos de los ejes y de la leyenda en sí según la Figura 4 del artículo.
plot(x = lions_merge$age, y = lions_merge$prop.black, frame = FALSE, main="Scatterplot", xlab = "Age (yr)", ylab = "Proportion black", pch = c(1, 16, 2, 17)[as.factor(lions_merge$group)])
legend('bottomright',legend=c("Serengeti females (n=62)","Serengeti males (n=22)","Ngorongoro females (n=11)","Ngorongoro males (n=10)"), pch=c(16,17,1,2))
(b) En el artículo se destacan los siguientes resultados: “After controlling for age, there was no effect of sex on nose colour in the Serengeti, but Ngorongoro males had lighter noses than Ngorongoro females”. Ajustar un primer modelo sin considerar la posible interacción entre el sexo y las áreas y contrastar si el sexo es significativo en el modelo así ajustado y en los modelos separados según el área.
En este caso, la variable dependiente sería la proporción del color de la nariz, puesto que es la que estaremos analizando en base a las demás. Se comenta que controlaron la edad como covariable, por lo que la tendremos en cuenta en los últimos modelos. Para este modelo inicial, como dicta el enunciado, no la incluiremos como variable independiente. Finalmente, debemos tener en cuenta que tanto el sexo como el área son variables categóricas dicotómicas, por lo que tendremos que saber cuál es tomada como “positiva” y cuál como “negativa” al analizar los valores estimadores. Para ello transformaremos las variables categóricas a factores y lo revisaremos con el comando “contrasts”.
lions$sex <- as.factor(lions$sex)
lions$area <- as.factor(lions$area)
contrasts(lions$sex)## M
## F 0
## M 1contrasts(lions$area)## S
## N 0
## S 1Vemos que se toma como 1 el sexo masculino y el área Serengeti.
Ahora generaremos primeramente un modelo de regresión lineal según la frase del enunciado, sin adjuntar la edad.
model_lions1 <- lm(prop.black ~ sex + area, data=lions)
summary(model_lions1)##
## Call:
## lm(formula = prop.black ~ sex + area, data = lions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48591 -0.19563 -0.02563 0.24409 0.47538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.52489 0.06430 8.163 9.03e-13 ***
## sexM -0.21028 0.05752 -3.655 0.000408 ***
## areaS 0.01101 0.06620 0.166 0.868208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2666 on 102 degrees of freedom
## Multiple R-squared: 0.1215, Adjusted R-squared: 0.1043
## F-statistic: 7.053 on 2 and 102 DF, p-value: 0.001352Para este primer modelo, el cual parece explicable (al tener un p-valor inferior al nivel de significación alfa = 0.05), observamos cómo el sexo es significativo (p-valor 0.000408), mientras que el área no lo es (0.868208, tomando la misma referencia de alfa). Como los predictores son variables categóricas (sex y area), utilizaremos el comando anova del la librería car (ya que automáticamente considera diseños no balanceados) para comprobar si estos resultados pueden ser verífidicos.
library(car)
anova(model_lions1)## Analysis of Variance Table
##
## Response: prop.black
## Df Sum Sq Mean Sq F value Pr(>F)
## sex 1 1.0005 1.00045 14.0781 0.0002916 ***
## area 1 0.0020 0.00197 0.0277 0.8682077
## Residuals 102 7.2486 0.07106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Observamos cómo, efectivamente, parece que el sexo es significativo (al tener un p-valor muy inferior al nivel de significación alfa = 0.05), por lo que puede parecer que éste tenga un comportamiento distintivo en la coloración negra de las narices de los leones.
Ahora analizaremos el comportamiento de la misma variable dependiente, pero dividiendo entre áreas. Para ello dividiremos el dataframe en dos según el area, y realizaremos el mismo modelo con estos datos. Como dicta el artículo, tomaremos la variable edad como variable concomitante.
# Para Ngorongoro:
lions_ngoro <- lions[lions$area == "N",]
model_ngoro <- lm(prop.black ~ sex*age, data=lions_ngoro)
summary(model_ngoro)##
## Call:
## lm(formula = prop.black ~ sex * age, data = lions_ngoro)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15771 -0.08862 -0.02669 0.06724 0.25969
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.32531 0.14975 -2.172 0.0443 *
## sexM 0.35005 0.19249 1.819 0.0866 .
## age 0.15848 0.03112 5.092 9.04e-05 ***
## sexM:age -0.10024 0.03498 -2.866 0.0107 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1293 on 17 degrees of freedom
## Multiple R-squared: 0.6991, Adjusted R-squared: 0.6461
## F-statistic: 13.17 on 3 and 17 DF, p-value: 0.000108# Para Serengeti
lions_seren <- lions[lions$area == "S",]
model_seren <- lm(prop.black ~ sex*age, data=lions_seren)
summary(model_seren)##
## Call:
## lm(formula = prop.black ~ sex * age, data = lions_seren)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30712 -0.08717 0.02071 0.09153 0.33028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.074893 0.035288 2.122 0.0369 *
## sexM -0.130055 0.076583 -1.698 0.0934 .
## age 0.075804 0.004905 15.454 <2e-16 ***
## sexM:age 0.031156 0.021058 1.480 0.1429
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1307 on 80 degrees of freedom
## Multiple R-squared: 0.8116, Adjusted R-squared: 0.8046
## F-statistic: 114.9 on 3 and 80 DF, p-value: < 2.2e-16Obtenemos dos modelos a priori válidos estadísticamente (con p-valores inferiores al nivel de significación estándar y con “Multiple r-squared” más próximas a 1 que a 0). No obstante, se muestran unos resultados que son distintos a los de la frase del enunciado. En este caso, al realizar un modelo únicamente con los leones de Ngorongoro y la interacción, encontramos que parece no existir una diferencia singificativa entre el sexo de los leones de la región (p-valor 0.0866). Por otro lado, tanto la edad como la interacción entre edad y sexo son significativas, por lo que parece que la edad de los leones de Ngorongoro puede afectar la coloración de su nariz, así como la interacción con su sexo. Para los leones en Serengeti, tampoco observamos diferencias significativas en el modelo individual, a nivel de sexo (p-valor 0.0934), pero sí observamos diferencias significativas en su edad.
Probaremos a analizar los modelos sin interacción entre áreas, para testear si la diferencia entre estos resultados y el enunciado se da por la interacción.
# Para Ngorongoro:
model_ngoro2 <- lm(prop.black ~ sex, data=lions_ngoro)
summary(model_ngoro2)##
## Call:
## lm(formula = prop.black ~ sex, data = lions_ngoro)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34000 -0.15091 -0.06091 0.08909 0.48909
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.41091 0.06709 6.125 6.9e-06 ***
## sexM 0.02909 0.09722 0.299 0.768
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2225 on 19 degrees of freedom
## Multiple R-squared: 0.00469, Adjusted R-squared: -0.04769
## F-statistic: 0.08954 on 1 and 19 DF, p-value: 0.768# Para Serengeti
model_seren2 <- lm(prop.black ~ sex, data=lions_seren)
summary(model_seren2)##
## Call:
## lm(formula = prop.black ~ sex, data = lions_seren)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50613 -0.16051 -0.04238 0.25137 0.45136
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.55613 0.03410 16.310 < 2e-16 ***
## sexM -0.28749 0.06663 -4.315 4.43e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2685 on 82 degrees of freedom
## Multiple R-squared: 0.185, Adjusted R-squared: 0.1751
## F-statistic: 18.62 on 1 and 82 DF, p-value: 4.428e-05Para los leones en Ngorongoro no obtenemos diferencias significativas a nivel de sexo (p-valor 0.768). Este modelo no es significativo (p-valor 0.768). No obstante, para los leones en Serengeti, parece que existe una diferencia significativa entre los leones macho y hembra de la región (con pvalor 4.43e-05), en un modelo que sí que parece ser significativo (p-valor 4.428e-05). Según este modelo, como el estimador por sexo es negativo (-0.28749) y el nivel de sexo positivo (1) es traducido como “macho”, parece que los machos tienen menos proporción de color negro en la nariz que las hembras en Serengeti. Aun así, el “multiple R-squared” del modelo es realmente bajo (mucho más cercano a 0 que 1), por lo que en realidad muy poca variabilidad de la proporción de negro en la nariz es explicada por el sexo. Como ambos modelos no son explicativos (al tener un R^2 muy próximas a cero) desestimaremos estos datos y tomaremos los datos anteriores como válidos.
(c) Otro resultado destacado es que para los machos hay diferencias según el área. Contrastar este resultado y dibujar las rectas de regresión para las dos áreas que se obtienen del modelo.
Para contrastar si existe diferencias dentro de los machos en las dos áreas, crearemos un dataframe con únicamente los machos. Como en el enunciado no se especifica si añadir la edad como variable independiente, probaremos un modelo con ella y otro sin ella, y decidiremos cuál es más explicativo.
Para empezar crearemos el modelo de regresión lineal con el porcentaje de negro en la nariz como variable dependiente, y el área como variable independiente.
lions_male <- lions[lions$sex == "M",]
model_male1 <- lm(prop.black ~ area, data=lions_male)
summary(model_male1)##
## Call:
## lm(formula = prop.black ~ area, data = lions_male)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34000 -0.12114 -0.05864 0.04750 0.45136
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.44000 0.05832 7.545 2.06e-08 ***
## areaS -0.17136 0.07033 -2.436 0.021 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1844 on 30 degrees of freedom
## Multiple R-squared: 0.1652, Adjusted R-squared: 0.1374
## F-statistic: 5.936 on 1 and 30 DF, p-value: 0.02098Con este primer modelo, podemos ver cómo el modelo anterior parece significativo (p-valor de 0.02098 respecto el nivel de significación alfa = 0.05). Éste anuncia que existe una diferencia significativa (con un p-valor de 0.021) entre las áreas estudiadas para los machos. Como su estimador es negativo (-0.17136) y teniendo en cuenta que el área tomada como positiva (1) es Serengeti, podemos deducir que en esta área los machos tienen la nariz más clara (tienen menos proporción de negro) que en Ngorongoro. No obstante, la R^2 del modelo es considerablemente baja (0.1652), mostrando que quizás el ajuste del modelo con estos datos no sea elevado.
Si añadimos la edad:
model_male2 <- lm(prop.black ~ age + area, data=lions_male)
summary(model_male2)##
## Call:
## lm(formula = prop.black ~ age + area, data = lions_male)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16192 -0.08356 -0.01158 0.08842 0.22278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.10826 0.08831 -1.226 0.2301
## age 0.07689 0.01126 6.827 1.7e-07 ***
## areaS 0.14411 0.06402 2.251 0.0321 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1162 on 29 degrees of freedom
## Multiple R-squared: 0.6798, Adjusted R-squared: 0.6577
## F-statistic: 30.78 on 2 and 29 DF, p-value: 6.745e-08Este segundo modelo nos da una R^2 mucho más elevada (0.6798), posibilemente por el alto nivel de explicación de la variable edad (la cual es significativa, con un p-valor de 1.7e-07). Como el modelo en sí es significativo (p-valor 6.745e-08), parece que éste es más explicativo y se ajusta más a los datos que el anterior. Pero como deseamos ver el posible efecto diferencial del área en los machos, la conclusión es la misma: sí que parece haber una diferencia significativa (en este caso p-valor de 0.0321) en la coloración de la nariz de los machos entre las dos áreas de estudio. La diferencia, no obstante, reside en que, para este modelo, el estimador es positivo, mientras que para el modelo anterior es negativo. Por lo tanto, existe un efecto potencial de la edad que hace modificar el output de este último modelo.
Para representar las rectas de regresión, dividiremos los datos según el área, aunque la representación se hará en conjunto, sobre los datos completos. Se tomará el modelo que incluye la edad como variable independiente. La recta para los machos de Ngorongoro se representará en azul, mientras que la recta para los machos de Serengeti se representará en rojo.
lions_male_ngoro <- lions_male[lions_male$area == "N",]
lions_male_seren <- lions_male[lions_male$area == "S",]
plot(prop.black ~ age, lions_male)
abline(recta1 <- lm(prop.black ~ age, lions_male_ngoro), col="blue")
abline(recta2 <- lm(prop.black ~ age, lions_male_seren), col="red")
(d) En la tabla 1 del artículo de Whitman et al. se dan los intervalos de confianza al 95 %, al 75% y al 50% para predecir la edad de una leona de 10 años o menos según su proporción de pigmentación oscura en la nariz. La primera cuestión es: ¿sirven para esto los modelos estudiados en los apartados anteriores?
El primer modelo (2b) tenía como variable dependiente la proporción de negro en la nariz, y el área y el sexo como variables independientes.
El segundo y tercer modelo realizados (2b) tenían como variable dependiente la proporción de negro en la nariz y como variable independiente, el sexo (dividido por áreas).
El cuarto modelo (2c) tenía como variable dependiente la proporción de negro en la nariz y como variable independiente, el área (teniendo en cuenta únicamente los machos).
Para esta tabla se ha utilizado la proporción de negro en la nariz para predecir la edad. Así, la proporción de negro en la nariz sería una de las variables independientes, mientras que la edad sería la variable dependiente. Por lo tanto, no son del todo útiles los modelos estudiados, ya que parten de premisas diferentes y, consecuentemente, responden a preguntas distintas.
Reproducir la fila de la tabla 1 para una proporción del 0.50 según el modelo que proponen en el artículo.
Según la descripción de la Tabla 1 del artículo, el modelo para este apartado se hizo con los mínimos cuadrados únicamente de las 63 hembras de ambas áreas que tengan 10 o menos años. Por lo tanto, filtraremos el dataframe original en base a las premisas que nos describen.
Para el modelo, ahora tendremos en cuenta que se quiere deducir la edad, por lo que será la variable dependiente. Como variables independientes anotaremos la proporción de negro en la nariz, puesto que el sexo es invariable en este caso (son todo hembras), y no se nos especifica de qué área proviene el individuo a predecir.
# Filtramos el dataframe
lions_female <- lions[lions$sex == "F",]
lions_female <- lions_female[lions_female$age <= 10,]
# Verificamos que la n es de 63, como dice la descripción de la Tabla 1.
nrow(lions_female)## [1] 63# Y creamos el modelo:
model_lionsage <- lm(age ~ prop.black, data=lions_female)
# Añadimos la nueva individuo (hembra, 0.5 proporción de negro)
newlion <- data.frame(prop.black=0.5, sex=F)
# Precedimos los tres intervalos de confianza
predict(model_lionsage, newlion, interval="prediction", level=0.95) ## fit lwr upr
## 1 5.381505 3.017659 7.745351predict(model_lionsage, newlion, interval="prediction", level=0.75) ## fit lwr upr
## 1 5.381505 4.00855 6.75446predict(model_lionsage, newlion, interval="prediction", level=0.5) ## fit lwr upr
## 1 5.381505 4.57938 6.183629Observamos que los valores que se han predecido són parecidos a los de la tabla (diferencias de unos 0.4 puntos, probablemente dadas por el uso de aproximados que se menciona en el enunciado).
Para este modelo tenemos que:
El intervalo de confianza para la edad de esta leona, al 95%, sería (3.018-7.75).
El intervalo de confianza para la edad de esta leona, al 75%, sería (4.009-6.75).
El intervalo de confianza para la edad de esta leona, al 95%, sería (4.58-6.18).
En cuanto al valor predecido de la edad, simplemente debemos sustiuir el valor de la proporción de negro (en este caso, 0.5) en la fórmula de la recta de regresión que proponen:
y = 2.0667 + 5.9037arcsin(0.5) = 2.0667 + 5.9037*0.5235987756 = 5.15787
Aclarar un detalle: lo que en la tabla 1 del artículo se llama s.e., standard error ¿qué es exactamente? Nota: Recordemos también aquí que los resultados pueden ser ligeramente distintos a los del artículo por la utilización de datos aproximados.
El error estándar que anotan consiste en el abarcamiento máximo y mínimo que puede abarcar la predicción, respecto la media calculada, sin tener en cuenta los intervalos de confianza.
Para calcular la predicción han utilizado la fórmula de la recta de regresión, la cual viene dada por un modelo de regresión lineal que, como la mayoría de modelos del estilo, tiene una variable de error (ε) intrínseco, que, al calcularla permite tener en cuenta la desviación respecto la media.
A modo de ejemplo, la primera predicción para una leona con un 10% de negro en la nariz (proporción 0.1) es 2.66 (con s.e., 1.24). Es decir, se espera que tenga 2.66 +/- 1.24 años (de 1.42 a 3.9 años).
Ejercicio 3
(a) Verificar las hipótesis de Gauss-Markov y la normalidad de los residuos del modelo completo del apartado (b) del ejercicio 2. Realizar una completa diagnosis del modelo para ver si se cumplen las condiciones del modelo de regresión: normalidad, homocedasticidad,. . . y estudiar la presencia de valores atípicos, de alto leverage y/o puntos influyentes. Construir los gráficos correspondientes y justificar su interpretación. ¿Podemos considerar el modelo ajustado como fiable?
De todos los modelos estudiados en el apartado b, intuimos que el modelo pedido es el siguiente, donde se toma la proporción de color negro de la nariz como variable dependiente y el sexo y el área como dependientes, puesto que los otros se crearon al dividir los sujetos por áreas.
Así pues, el modelo inicial será:
model_lions1 <- lm(prop.black ~ sex + area, data=lions)Las condiciones de Gauss-Markov son:
- Que los errores tienen a una media teórica (esperanza) de cero. Es decir:
E(ϵi) = 0 i = 1, . . . , n
- La varianza σ2 de todos los errores es la misma. Es decir:
var(ϵi) = E(ϵ2i) = σ2 i = 1, . . . , n
- Los errores son incorrelacionados entre sí. Es decir:
E(ϵiϵj) = 0 para todo i ̸= j
Por lo tanto, necesitaremos contrastar la información de los residuos. Estos los podemos observar mediante el comando “residuals”, comprobar que son muy próximos a cero, y que su media, por lo tanto, también lo es. Mediante este cálculo podemos intuir si existe linealidad en los datos.
Representaremos los residuos en un histograma (para visualizarlos a nivel general) y en un boxplot (para enfatizar la media de los mismos y su distribución respecto el cero).
# Cálculo de residuos y su media:
residuals(model_lions1)## 1 2 3 4 5 6
## -0.115628761 -0.185628761 -0.215628761 -0.195628761 -0.205628761 -0.195628761
## 7 8 9 10 11 12
## -0.205628761 -0.145628761 -0.095628761 -0.105628761 -0.125628761 -0.155628761
## 13 14 15 16 17 18
## -0.175628761 -0.055628761 -0.065628761 -0.115628761 -0.025628761 0.094371239
## 19 20 21 22 23 24
## 0.104371239 0.264371239 0.274371239 0.394371239 -0.024616725 -0.214616725
## 25 26 27 28 29 30
## 0.165383275 0.125383275 0.025383275 0.055383275 0.025383275 0.425383275
## 31 32 33 34 35 36
## 0.475383275 0.195383275 -0.485905923 -0.455905923 -0.475905923 -0.385905923
## 37 38 39 40 41 42
## -0.395905923 -0.475905923 -0.345905923 -0.455905923 -0.385905923 -0.325905923
## 43 44 45 46 47 48
## -0.315905923 -0.285905923 -0.245905923 -0.165905923 -0.175905923 -0.295905923
## 49 50 51 52 53 54
## -0.295905923 -0.215905923 -0.185905923 -0.145905923 -0.105905923 -0.085905923
## 55 56 57 58 59 60
## -0.045905923 -0.025905923 -0.005905923 -0.055905923 -0.175905923 -0.135905923
## 61 62 63 64 65 66
## -0.075905923 -0.005905923 0.024094077 -0.055905923 0.094094077 0.054094077
## 67 68 69 70 71 72
## 0.084094077 0.154094077 0.164094077 0.244094077 0.324094077 0.154094077
## 73 74 75 76 77 78
## 0.094094077 0.224094077 0.314094077 0.314094077 0.294094077 0.264094077
## 79 80 81 82 83 84
## 0.264094077 0.294094077 0.324094077 0.354094077 0.374094077 0.364094077
## 85 86 87 88 89 90
## 0.454094077 0.454094077 0.414094077 0.444094077 0.324094077 0.294094077
## 91 92 93 94 95 96
## 0.424094077 0.454094077 0.384094077 0.124094077 -0.364893887 -0.374893887
## 97 98 99 100 101 102
## -0.304893887 -0.264893887 -0.274893887 -0.024893887 0.075106113 0.105106113
## 103 104 105
## 0.375106113 -0.024893887 -0.174893887mean(model_lions1$residuals)## [1] -2.27463e-17stem(residuals(model_lions1))##
## The decimal point is 1 digit(s) to the left of the |
##
## -4 | 98866
## -4 | 0
## -3 | 99765
## -3 | 32000
## -2 | 9765
## -2 | 2211100
## -1 | 9988877655
## -1 | 4322110
## -0 | 9876665
## -0 | 3322211
## 0 | 233
## 0 | 5688999
## 1 | 0123
## 1 | 5567
## 2 | 024
## 2 | 6667999
## 3 | 11222
## 3 | 567889
## 4 | 1234
## 4 | 5558# Los representamos:
hist(residuals(model_lions1),xlab="Residuos", ylab="Frecuencia", main="Residuos del modelo 2b (model_lions1)")
boxplot(residuals(model_lions1), main="Residuos del modelo 2b (model_lions1)")
En el histograma observamos que los residuos no siguen una distribución normal perfecta, habiendo un número mayor de errores negativos que desplazan el pico de la normal hacia la izquierda. Por otro lado, en el boxplot se nos informa de que la mayor parte de errores tiene valores cercanos a cero. Más adelante comprobaremos si esta hipótesis sobre la no-normalidad se cumple.
Para determinar si existe homocedasticidad o heterocedasticidad de los datos, podemos utilizar un test como el de Breusch-Pagan. Para éste, la hipótesis nula (Ho) será que existe homocedasticidad en los residuos, mientras que la hipótesis alternativa (H1) será que existe heterocedasticidad para ellos.
library(lmtest)## Warning: package 'lmtest' was built under R version 4.2.2## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.2.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(model_lions1)##
## studentized Breusch-Pagan test
##
## data: model_lions1
## BP = 8.8957, df = 2, p-value = 0.0117Con un pvalor inferior al nivel de significación estándar (alfa=0.05) rechazamos la hipótesis nula y aceptamos que existe heterocedasticidad en los datos (es decir, no se cumple la homocedasticidad). Por lo tanto, parece que la varianza de los errores no es igual en todas las observaciones realizadas.
Continuaremos el diagnóstico mediante gráficos que nos ayuden a visualizar el comportamiento de los residuos:
Gráfico de dispersión de los residuos respecto al índice
Gráfico de residuos vs valores ajustados
Gráfico de valores observados vs valores ajustados
Gráfico QQ-plot
par(mfrow=c(2,2))
plot(residuals(model_lions1),xlab="Índice",ylab="Residuos", main = "Dispersión de Residuos vs Índice")
plot(fitted(model_lions1),residuals(model_lions1),xlab="Ajustados",ylab="Residuos", main ="Residuos vs. Valores ajustados")
plot(fitted(model_lions1),lions$prop.black,xlab="Valores ajustados",ylab="Valores observados", main="Valores observados vs valores ajustados")
abline(0,1)
qqnorm(residuals(model_lions1),xlab="Cuantiles de distribución normal",ylab="Residuos", main="QQ-plot")
qqline(residuals(model_lions1))
par(mfrow=c(1,1))El gráfico de dispersión de los residuos respecto al índice muestra si existe alguna agrupación que difiere de la aleatoriedad. En este caso debemos tener en cuenta que en el modelo tenemos únicamente variables categóricas dicotómicas como variables explicativas (sexo y área) y de ahí que pueda observarse el patrón diagonal de los residuos respecto el índice.
Al tratar con variables categóricas, en el gráfico de residuos vs valores ajustados observamos ambas variables independientes, con sus dos niveles cada una (cuatro líneas verticales). Debemos fijarnos en este caso que los residuos (círculos del gráfico) estén dispersos a lo largo de la línea vertical, de forma que se cumpla la aleatoriedad. En este caso parece cumplirse aeatoriedad excepto por el primer y segundo grupo, donde se tienden a acumular los residuos en vez de estar repartidos.
El gráfico de valores observados vs valores ajustados muestra el ajuste de la recta de regresión. Para ello, debemos observar la proximidad de los puntos a la bisectriz (línea recta dibujada). Vemos cómo la línea no es horizontal, reafirmando que los datos no siguen un comportamiento aleatorio.
Como tratamos con variables categóricas, para este paso consideramos más adecuado estudiar el ajuste del modelo de otra forma, como por ejemplo mediante el test de ajuste que nos muestra el mismo “summary” del modelo, comentado anteriormente.
- Para observar la normalidad de los residuos, podemos analizar su QQ-plot. En éste, los círculos representan los puntos asociados a los cuantiles de una distribución normal (por defecto, estándar), mientras que la recta representa el ajuste perfecto a esta distribución normal. Como podemos observar, los círculos siguen el comportamiento de la recta en su mayoría, con algunos valores dispersos en los extremos. Como no queda claro, sería adecuado utilizar un test de normalidad, como el de Shapiro-Wilk, para verificarlo. La hipótesis nula (Ho) de este test es que existe normalidad de los datos (en este caso, de los residuos). La hipótesis alternativa (H1) es que NO existe normalidad en los datos (es decir,éstos no siguen una distribución normal).
shapiro.test(residuals(model_lions1))##
## Shapiro-Wilk normality test
##
## data: residuals(model_lions1)
## W = 0.96366, p-value = 0.005624Con un p-valor inferior al nivel de significación estándar alfa = 0.05, rechazamos la hipótesis nula, por lo que, según el test, parece que NO existe normalidad en los residuos. Estos datos son coherentes con la distribución que hemos observado en el primer plot.
El cálculo del leverage lo realizaremos mediante los valores de leverage (el comando hatpoints). Los ordenaremos para visibilizar aquellos más elevados y los representaremos en un gráfico de barras.
leverage_points <- as.data.frame(hatvalues(model_lions1))
leverage_points[order(-leverage_points['hatvalues(model_lions1)']), ]## Warning in xtfrm.data.frame(x): cannot xtfrm data frames## [1] 0.06039489 0.06039489 0.06039489 0.06039489 0.06039489 0.06039489
## [7] 0.06039489 0.06039489 0.06039489 0.06039489 0.05817759 0.05817759
## [13] 0.05817759 0.05817759 0.05817759 0.05817759 0.05817759 0.05817759
## [19] 0.05817759 0.05817759 0.05817759 0.03727167 0.03727167 0.03727167
## [25] 0.03727167 0.03727167 0.03727167 0.03727167 0.03727167 0.03727167
## [31] 0.03727167 0.03727167 0.03727167 0.03727167 0.03727167 0.03727167
## [37] 0.03727167 0.03727167 0.03727167 0.03727167 0.03727167 0.03727167
## [43] 0.03727167 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [49] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [55] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [61] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [67] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [73] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [79] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [85] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [91] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [97] 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872 0.01509872
## [103] 0.01509872 0.01509872 0.01509872barplot(hatvalues(model_lions1))
# Podemos cuantificar cuántos están por encima de 0.05, como guía, y calcular su porcentaje.
length(which(hatvalues(model_lions1)>0.05))/length(hatvalues(model_lions1))*100## [1] 20# Anotamos los índices de los valores en cuestión
which(hatvalues(model_lions1)>0.05)## 23 24 25 26 27 28 29 30 31 32 95 96 97 98 99 100 101 102 103 104
## 23 24 25 26 27 28 29 30 31 32 95 96 97 98 99 100 101 102 103 104
## 105
## 105Existe un 20% de puntos de alto leverage, tomando el valor de 0.05 como referencia. Estos se desvían del comportamiento esperado y pueden ser puntos influyentes potenciales.
Seguiremos analizando posibles puntos influyentes mediante el cálculo de la distancia de Cook (Cook y Weisberg 1994) para cada valor anotado del modelo. Luego, realizaremos la representación con un gráfico de dispersión. Según la bibliografía, para una n de 105, una distancia de Cook cuatro veces superior a la distancia media implica la presencia de puntos influyentes. Por lo tanto, representaremos la línea que separa esta distancia, y enumeraremos aquellos puntos que sean influyentes.
library(car)
# Cálculo de la distancia de Cook.
cooksd <- cooks.distance(model_lions1)
# Representación
plot(cooksd, pch="*", cex=2, main="Puntos influyentes (Cook)")
# Determinación de la línea umbral
abline(h = 4*mean(cooksd, na.rm=T), col="red")
# Anotamos los puntos con distancia de Cook superior a 4 veces la media de todos los puntos
text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),""), col="red")
Observamos que el índice de todos estos puntos coincide con el índice de algunos de los veinte puntos con alto leverage anotados anteriormente (puntos 30, 31, 95, 96 y 103). Por lo tanto, parece lógico que estos valores son outliers y que, probablemente, deberían ser eliminados para una optimización del modelo.
Finalmente, para responder sobre la fiabilidad de este modelo, añadiremos el apunte de su R^2.
summary(model_lions1)##
## Call:
## lm(formula = prop.black ~ sex + area, data = lions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48591 -0.19563 -0.02563 0.24409 0.47538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.52489 0.06430 8.163 9.03e-13 ***
## sexM -0.21028 0.05752 -3.655 0.000408 ***
## areaS 0.01101 0.06620 0.166 0.868208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2666 on 102 degrees of freedom
## Multiple R-squared: 0.1215, Adjusted R-squared: 0.1043
## F-statistic: 7.053 on 2 and 102 DF, p-value: 0.001352Con un R^2 de 0.1215, realmente lejana al 1 utópico, unos residuos que se alejan de una distribución normal y diversos outliers y puntos influyentes, podemos afirmar que este modelo no es del todo fiable.
(b) Teniendo en cuenta que la variable respuesta de la regresión del apartado (b) del ejercicio 2 es una proporción, ¿presenta algún problema este modelo? ¿Qué alternativas nos podemos plantear para mejorar el ajuste de los datos? Atendiendo a la naturaleza de la variable respuesta, ¿hay alguna transformación adecuada?
El hecho que sea una proporción, indica que la propia variable tiene un límite determinado (es decir, va de 0 a 1). Esto acota mucho la distribución de los datos, y hace que el ajuste del modelo pueda verse involucrado (puede ser una de las razones por la cual la R^2 es tan baja).
Para optimizar el ajuste de los datos podemos realizar una transformación de los mismos, además de quitar los puntos influyentes (outliers, puntos de alto leverage…) que hemos detectado anteriormente.
Optando por la opción de la transformación, probaremos de realizar una transformación arcoseno de la proporción de color negro en la nariz, que, según la webgrafía, es la adecuada para este tipo de variables que abarcan de 0 a 1, o bien porcentajes.
(c) Aplicar la transformación más adecuada a la variable respuesta del modelo considerado. Comparar los dos modelos: con y sin la transformación. ¿Qué modelo es mejor? Justificar la respuesta.
Para la transformación, utilizaremos el arcoseno de la variable en cuestión (prop.black) y las variables independientes serán las mismas. Mediante el summary veremos los estadísticos base del modelo y lo compararemos con el summary del modelo anterior.
model_lions_asin <- lm(asin(prop.black) ~ sex + area, data=lions)
summary(model_lions_asin)##
## Call:
## lm(formula = asin(prop.black) ~ sex + area, data = lions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57444 -0.23034 -0.06586 0.27020 0.80480
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.58791 0.08478 6.935 3.8e-10 ***
## sexM -0.27451 0.07585 -3.619 0.000462 ***
## areaS 0.03655 0.08728 0.419 0.676250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3515 on 102 degrees of freedom
## Multiple R-squared: 0.1232, Adjusted R-squared: 0.106
## F-statistic: 7.168 on 2 and 102 DF, p-value: 0.001222summary(model_lions1)##
## Call:
## lm(formula = prop.black ~ sex + area, data = lions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48591 -0.19563 -0.02563 0.24409 0.47538
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.52489 0.06430 8.163 9.03e-13 ***
## sexM -0.21028 0.05752 -3.655 0.000408 ***
## areaS 0.01101 0.06620 0.166 0.868208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2666 on 102 degrees of freedom
## Multiple R-squared: 0.1215, Adjusted R-squared: 0.1043
## F-statistic: 7.053 on 2 and 102 DF, p-value: 0.001352Si nos fijamos en el p-valor del modelo con escala arcosénica, observamos que los resultados del modelo parecen estadísticamente significativos (p-valor de 0.001222, inferior al nivel de significación alfa = 0.05). Esta caraterística también se daba en el modelo inicial (p-valor de 0.001352).
Para comparar la bondad de ajuste, una vez más nos fijaremos en la R^2 de ambos modelos. Para el modelo con escala arcosénica, éste tiene una multiple R-squared de 0.1232, mientras que para el modelo inicial era de 0.1215. Este aumento ligero en la R^2 puede implicar que el modelo en escala arconsénica está mejor ajustado que el modelo inicial, aunque en ambos casos siguen siendo muy distantes a 1, intuyendo que se podría probar en realizar algún otro tipo de transformación, probablemente.
En cuanto a los resultados, la información que nos otorgan ambos modelos es muy parecida. Ambos determinan que existen diferencias significativas en la proporción de color negro de la nariz, o bien el arcoseno de la misma, según el sexo. En el caso del modelo en escala arcosénica, su estimador es ligeramente mayor para esta variable (-0.27451). Aun así, viendo el valor tan pequeño de la R^2, ponemos en duda su veracidad, al faltarle un buen ajuste.
(d) Realizar una rápida diagnosis del modelo transformado. ¿Estamos satisfechos con este nuevo modelo? ¿Qué otro ajuste nos podemos plantear para mejorar el modelo?
Como antes, visualizaremos el comportamiento de los residuos.
plot(residuals(model_lions_asin), ylab="Residuos")
abline(h=0)
Observamos como el arcoseno no ha cambiado demasiado la distribución de los residuos. Anteriormente parecían existir más residuos con valor positivo, de forma que se ha compensado ligeramente la afectación de este hecho. No obstante, el comportamiento de los residuos dista del óptimo, alejándose de cero.
Para analizar la normalidad, homocedasticidad y otras características de los datos en escala arcoseno, realizaremos diferentes plots del modelo.
plot(model_lions_asin)
En el gráfico de los residuos vs valores ajustados observamos, como antes, ambas variables independientes, con sus dos niveles cada una (cuatro líneas verticales). Los residuos (círculos del gráfico) siguen dispersos a lo largo de su línea vertical por lo que existe cierta aleatoriedad. No obstante, siguen habiendo zonas de acumulación de residuos, haciendo que la recta roja de comportamiento general sea distante a una línea horizontal. En el gráfico 3, de escalado, se realiza la raíz cuadrada de los residuos estandarizados (eje y), permitiendo que el comportamiento de los mismos sea ligeramente más aleatorio, aunque la recta resultante dista de ser horizontal.
El gráfico de residuos vs leverage muestra los puntos de alto leverage y los enumera con su índice. Como antes, para ello se ha calculado la distancia de Cook, y observamos que existe una cierta desviación de la línea de guía (línea roja en comparación con la discontinua), debido a estos puntos. En este caso solamente se marcan tres puntos (30, 31 y 103), que están incluidos en la lista que se ha realizado en el ejercicio anterior con puntos de alto leverage, por lo que parece coherente.
Finalmente, en cuanto al QQplot de normalidad, observamos cómo los residuos estandarizados en el modelo arcosénico siguen la línea de normalidad, alejándose en los extremos. Este plot nos enumera aquellos posibles outliers en la parte superior derecha (85, 86 y 92). Como antes, como no queda claro visualmente, podemos utilizar un test estadístico para evaluar la normalidad.
shapiro.test(residuals(model_lions_asin))##
## Shapiro-Wilk normality test
##
## data: residuals(model_lions_asin)
## W = 0.96157, p-value = 0.003919Con un p-valor inferior al nivel de significación estándar alfa = 0.05 (0.003919), rechazamos la hipótesis nula, por lo que, según el test, no existe normalidad en los residuos.
Por lo tanto, parece que el modelo arcosénico dista de ser el óptimo.
Otra alternativa sería la transformación logarítmica, aunque se aplica normalmente para estabilizar varianzas. Probémoslo.
model_lions_log <- lm(log(prop.black) ~ sex + area, data=lions)
summary(model_lions_log)##
## Call:
## lm(formula = log(prop.black) ~ sex + area, data = lions)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1378 -0.5283 0.1240 0.6095 1.0145
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.76788 0.17379 -4.418 2.48e-05 ***
## sexM -0.48236 0.15548 -3.102 0.00248 **
## areaS -0.09007 0.17892 -0.503 0.61576
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7205 on 102 degrees of freedom
## Multiple R-squared: 0.08627, Adjusted R-squared: 0.06835
## F-statistic: 4.815 on 2 and 102 DF, p-value: 0.01004Aunque siga siendo significativo, para este modelo observamos que la R^2 decrece aún más, por lo que lo desestimamos directamente.
(e) Discutir la utilización de la transformación arcoseno en el modelo del apartado (d) del ejercicio 2.
El modelo mencionado en el apartado sigue la recta de regresión:
y = 2.0667 + 5.9037arcsin(x)
Donde X sigue siendo la proporción de coloración negra en la nariz de los leones. Como se ha comentado en el apartado anterior, una proporción es un valor limitado de 0 a 1, basado en un numerador de interés respecto el total. Por lo tanto, el abarcamiento de la variable es muy distinto respecto otras variables numéricas aleatorias, y el arcoseno permite corregir este tipo de desajuste.
Por lo tanto, parece una buena elección de transformación.
Webgrafía:
R. Dennis Cook, Sanford Weisberg (1994) An Introduction to Regression Graphics. John Wiley & Sons, Inc. DOI:10.1002/9780470316863
Regresión lineal
Gráficos
https://www.datanovia.com/en/blog/pch-in-r-best-tips/
http://www.sthda.com/english/wiki/r-plot-pch-symbols-the-different-point-shapes-available-in-r