Ejercicios modelo lineal. III

3. Ejercicios de Regresión Lineal Simple (segunda parte)

Ejercicio 3.1 Medidas del cuerpo, Parte IV. Base de datos bdims del paquete openintro.

1. Compare los ajustes realizados en los ejercicios 2.1 y 2.2.

• modelo 2.1:
  • formula = wgt ~ hip.gi
  • \(R^2\) 0.5821
• modelo 2.2:
  • formula = wgt ~ hgt
  • \(R^2\) 0.5145

El modelo wgt ~ hip.gi es mejor, porque tiene un mejor \(R^2\)

b)Para el ajuste del peso usando la circunferencia de cadera como única covariable, halle un intervalo de confianza de nivel 0.95 cuando el contorno de cadera mide 100 cm. Compárelo con el intervalo de predicción para ese mismo contorno de cadera.

lm_fitted <- lm(wgt~hip.gi,data =bdims )
new <- data.frame(hip.gi = 100)
#intervalo de confianza
predict.lm(lm_fitted,newdata = new,interval="confidence",level = 0.95)

       fit      lwr    upr
1 74.20646 73.36492 75.048

#intervalo de predicción
predict.lm(lm_fitted,newdata = new,interval="prediction",level = 0.95)

       fit      lwr      upr
1 74.20646 57.21927 91.19365

3. Para el ajuste del peso usando la altura como única covariable, halle un intervalo de confianza de nivel 0.95 cuando la altura es de 176 cm. Compárelo con el intervalo de predicción para esa misma altura. ¿Cuál de los dos modelos da un intervalo de predicción más útil?

lm_fitted <- lm(wgt~hgt,data =bdims )
new <- data.frame(hgt = 176)
#intervalo de confianza
predict.lm(lm_fitted,newdata = new,interval="confidence",level = 0.95)

      fit      lwr      upr
1 74.0893 73.17511 75.00349

#intervalo de predicción
predict.lm(lm_fitted,newdata = new,interval="prediction",level = 0.95)

      fit      lwr      upr
1 74.0893 55.77921 92.39939

Es más útil el primer modelo, de la circunferencia de la cadera, porque el intervalo de predicción, a un mismo nivel de significatividad, es más acotado.

4. Construya un intervalo de confianza para el peso esperado cuando el contorno de cintura es de 80cm.,95cm., 125cm. de nivel 0.95. Estos tres intervalos, ¿tienen nivel simultáneo 0.95? Es decir, la siguiente afirmación ¿es verdadera o falsa? Justifique. En aproximadamente 95 de cada 100 veces que yo construya los IC basados en una (misma) muestra, cada uno de los 3 IC contendrán al verdadero valor esperado del peso.

lm_fitted <- lm(wgt~hip.gi,data =bdims )
new <- data.frame(hip.gi = c(80,95,125))
#intervalo de confianza
predict.lm(lm_fitted,newdata = new,interval="confidence",level = 0.95)

        fit       lwr       upr
1  43.72306  41.69463  45.75148
2  66.58561  65.80858  67.36264
3 112.31071 109.02587 115.59554

Es falso, no tienen un nivel de confianza simultaneo del 95%, porque no son independientes los test, al tratarse de una misma muestra.

Ejercicio 3.2: heights.txt

1. Realice un scatterplot de los datos, con la altura de las madres en el eje horizontal.

library(alr3)

Loading required package: car

Attaching package: ‘car’

The following object is masked from ‘package:openintro’:

    densityPlot

The following object is masked from ‘package:dplyr’:

    recode

The following object is masked from ‘package:purrr’:

    some

height.plot <-  ggplot(heights, aes(Mheight,Dheight)) +
  geom_point()+
  geom_abline(colour="firebrick", size =1)+
  theme_minimal()+
  lims(x = c(55,75), y = c(55,75))
height.plot

¿Describe esta figura un buen resumen de la relación entre ambas variables? Si, aunque pareciera que existe una tendencia a que las hijas sean más altas que las madres.

En base al scatterplot, ¿parecería ser cierto que las madres más altas suelen tener hijas más altas y viceversa con las más bajas? Si b) Ajuste el modelo lineal a los datos. Indique el valor de la recta ajustada. Superpóngala al scatter plot.

lm.ajustado <- lm(Dheight~Mheight,data = heights)
lm.ajustado

Call:
lm(formula = Dheight ~ Mheight, data = heights)

Coefficients:
(Intercept)      Mheight  
    29.9174       0.5417  

#intervalo de confianza para la pendiente
height.plot+  geom_smooth(method = "lm", se = FALSE)

¿Presenta visualmente un mejor ajuste que la recta identidad postulada en el ítem anterior? si

#coeficientes, errores estandar, coeficientes de determinación
summary(lm.ajustado)

Call:
lm(formula = Dheight ~ Mheight, data = heights)

Residuals:
   Min     1Q Median     3Q    Max 
-7.397 -1.529  0.036  1.492  9.053 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 29.91744    1.62247   18.44   <2e-16 ***
Mheight      0.54175    0.02596   20.87   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.266 on 1373 degrees of freedom
Multiple R-squared:  0.2408,    Adjusted R-squared:  0.2402 
F-statistic: 435.5 on 1 and 1373 DF,  p-value: < 2.2e-16

# varianza de los errores
anova(lm.ajustado)

Analysis of Variance Table

Response: Dheight
            Df Sum Sq Mean Sq F value    Pr(>F)    
Mheight      1 2236.7 2236.66  435.47 < 2.2e-16 ***
Residuals 1373 7052.0    5.14                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

confint(lm.ajustado,level = 0.95,parm = 2)

            2.5 %    97.5 %
Mheight 0.4908201 0.5926739

El p-valor del test t para \(\beta_0\) es <2e-16, por lo que rechazamos la hipótesis nula de que \(E(Dheight |Mheight) = \beta_0\), esto implica que el incremento de la altura de la madre es significativa para la altura esperada de la hija

3. Prediga y obtenga un intervalo de predicción para la altura de una hija cuya madre mide 64 pulgadas.

new <- data.frame(Mheight = 64, Dheight = NA)
predict.lm(lm.ajustado,newdata = new,interval="confidence",level = 0.95)

       fit      lwr      upr
1 64.58925 64.44578 64.73271

4. Una pulgada equivale a 2.54cm. Convierta ambas variables a centímetros (Dheightcm y Mheightcm) y ajuste un modelo lineal a estas nuevas variables.

heights <- heights %>% 
  mutate(Mheightcm = Mheight/2.54,
         Dheightcm = Dheight/2.54)
lm.ajustado.cm <- lm(Dheightcm~Mheightcm,data = heights)
# Ajuste en cm
summary(lm.ajustado.cm)

Call:
lm(formula = Dheightcm ~ Mheightcm, data = heights)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9124 -0.6018  0.0142  0.5874  3.5640 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 11.77852    0.63877   18.44   <2e-16 ***
Mheightcm    0.54175    0.02596   20.87   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8922 on 1373 degrees of freedom
Multiple R-squared:  0.2408,    Adjusted R-squared:  0.2402 
F-statistic: 435.5 on 1 and 1373 DF,  p-value: < 2.2e-16

# Ajuste en pulgadas
summary(lm.ajustado)

Call:
lm(formula = Dheight ~ Mheight, data = heights)

Residuals:
   Min     1Q Median     3Q    Max 
-7.397 -1.529  0.036  1.492  9.053 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 29.91744    1.62247   18.44   <2e-16 ***
Mheight      0.54175    0.02596   20.87   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.266 on 1373 degrees of freedom
Multiple R-squared:  0.2408,    Adjusted R-squared:  0.2402 
F-statistic: 435.5 on 1 and 1373 DF,  p-value: < 2.2e-16

preguntas

• ¿Deberían cambiar los estimadores de \(\beta_0\) y \(\beta_1\)?
  • Si la ordenada al origen, \(\beta_0\), no \(\beta_1\) porque es independiente de las unidades de medida
• ¿De qué manera?
  • \(\beta_0\) tiene que ser en cm 2,54 veces menor que en pulgadas
• ¿Y los errores estándares?
  • de igual manera
• ¿Y los p-valores?
  • no cambian
• ¿Y el coeficiente de determinación?
  • no cambia
• ¿Y la estimación del desvío estándar de los errores?
  • tiene que ser en cm 2,54 veces menor que en pulgadas

Ejercicio 3.3. Simulación 1.

1. Generar n = 22 datos que sigan el modelo lineal simple \[ Y = 10 + 5X + \epsilon \]

\(\epsilon \sim \mathcal{N}(0,\,\sigma^{2}).\) con \(\sigma^2=49\)

set.seed(54321)
epsilon <- rnorm(n = 22,sd = 7)
densityPlot(epsilon, col="lightcoral")

X <- round(runif(n = 22, min = 0, max = 10),2)
Y <- 10 + 5*X + epsilon
XY <- data.frame(X=X, Y=Y)
ggplot(XY, aes(X,Y)) + geom_point()

3. Ajuste el modelo lineal

ajuste <- lm(Y~X, data = XY)
summary(ajuste)

Call:
lm(formula = Y ~ X, data = XY)

Residuals:
    Min      1Q  Median      3Q     Max 
-18.645  -7.530   1.272   7.506  16.880 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  13.2363     3.9821   3.324  0.00338 ** 
X             4.6735     0.7081   6.600 1.99e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.299 on 20 degrees of freedom
Multiple R-squared:  0.6853,    Adjusted R-squared:  0.6696 
F-statistic: 43.56 on 1 and 20 DF,  p-value: 1.986e-06

confint(ajuste)

               2.5 %    97.5 %
(Intercept) 4.929862 21.542819
X           3.196389  6.150661

preguntas:

• ¿Los verdaderos \(\beta_0\) y \(\beta_1\) pertenecen a dichos intervalos?
  • Si
• ¿Cuánto dio la pendiente estimada, \(\hat{\beta_1}\)?
  • 8.427
• ¿En qué parte de la salida del ajuste lineal podemos encontrar el estimador de \(\sigma\)? * en el Residual standard error se indica \(\sigma\)
• ¿Cuánto debería valer?
  • 7

4. Pídamosle al R que chequee si el 5 pertenece al IC de nivel 0.95 calculado en base a la muestra.

beta1est <- 5 > confint(ajuste, parm = 2)[1] & 5 < confint(ajuste, parm = 2)[2]
beta1est

[1] TRUE

5. Superpóngale al scatterplot de los datos la recta verdadera (en azul) y la estimada en base a ellos (en rojo).

ggplot(XY, aes(X,Y)) + 
  geom_point()+
  geom_smooth(method = "lm",color="firebrick",se = FALSE)+
  geom_abline(slope =5,intercept = 10, color = "steelblue", size = 1)

Ejercicio 3.4 Simulación 2.

1000 repeticiones del ejercicio anterior

1. ¿Puede usted anticipar, desde la teoría las respuestas de las preguntas que siguen?

1. Las pendientes estimadas en las B = 1000 replicaciones, ¿serán siempre iguales o cambiarán de replicación en replicación?

• Cambian

2. ¿Alrededor de qué número variarán las pendientes estimadas en las 1000 replicaciones?

• 5

3. Si hacemos un histograma de estas B = 1000 replicaciones, ¿a qué distribución debería parecerse?

• \(\beta_1 \sim \mathcal{N}(5,\,\sigma_b^{2}).\)

4. Aproximadamente, ¿qué porcentaje de los 1000 intervalos de confianza para la pendiente estimados a partir de las 1000 muestras cubrirá al verdadero valor de la pendiente?

• 950

beta1est <- rep(NA,1000)
icbet <- rep(NA,1000)
set.seed(54321)
X <- round(runif(n = 22, min = 0, max = 10),2)
for(i in 1:1000){
epsilon <- rnorm(n = 22,sd = 7)
Y <- 10 + 5*X + epsilon
XY <- data.frame(X=X, Y=Y)
ajuste <- lm(Y~X, data = XY)
beta1est[i] <- ajuste$coefficients[2]
icbet[i] <- 5 > confint(ajuste, parm = 2)[1] & 5 < confint(ajuste, parm = 2)[2]
}

4. Haga un histograma de las pendientes estimadas. ¿Qué distribución parecen tener los datos?

ggplot(as.data.frame(beta1est), aes(beta1est))+
  geom_histogram(fill='aquamarine3',color="gray40", alpha=0.75)+
  geom_vline(xintercept = mean(beta1est), linetype ="dashed")+
  theme_light()

Tienen una distribución aproximadamente normal, con media 5.0062792

5. ¿Qué proporción de los intervalos de conanza construidos contiene al verdadero valor de la pendiente?

sum(icbet)/length(icbet)

[1] 0.942

Ejercicio 3.5 Mamíferos, Parte IV.conjunto de datos mammals del paquete openintro.

1. Si dos animales difieren en el BodyWt por un factor de diez, dé un intervalo del 95% de confianza para la diferencia en el log10(BrainWt) para estos dos animales.

lm_fitted <- lm(logBrainWt ~ logBodyWt, data = mammals)
confint(lm_fitted)[2,]

    2.5 %    97.5 % 
0.6947503 0.8086215 

2. Para un mamífero que no está en la base de datos, cuyo peso corporal es de 100 kg., obtenga la predicción y un intervalo de nivel 95% de predicción del log10 (BrainWt) :

nuevo <- data.frame(logBodyWt = log10(100), logBrainWt = NA)
predict.lm(lm_fitted,newdata = nuevo,interval="prediction",level = 0.95)

       fit      lwr     upr
1 3.638161 2.237701 5.03862

intervalo de predicción del peso del cerebro del mamífero cuyo peso corporal es 100kg. Mirando los valores numéricos obtenidos, ¿parece muy útil el resultado obtenido?

10^predict.lm(lm_fitted,newdata = nuevo,interval="prediction",level = 0.95)[2:3]

[1]    172.8626 109299.9493

No es útil el resultado obtenido.

Ejercicio 3.6. wblake

1. Hacer un scatter plot de la longitud (Length) en función de la edad (Age).

ggplot(wblake, aes(Age, Length))+geom_point()

La variable años es discreta

2. Calcule las medias y los desvíos estándares muestrales para cada uno de las ocho subpoblaciones

wblake %>% 
  group_by(Age) %>% 
  summarise(media   = mean(Length),
            desvio_estandar = sd(Length))

ggplot(wblake, aes(Age, Length, fill = as.factor(Age), group = Age))+
  geom_boxplot()+
  theme_minimal()+
  theme(legend.position = "none")

La longitud, ¿parece aumentar con la edad? si
La dispersión de la longitud, ¿parece mantenerse más o menos constante con la edad? ? ¿O crece? ¿O decrece? se mantiene constante

3. Ajuste un modelo lineal para explicar la longitud (Length) en función de la edad (Age).

ajuste <- lm(Length ~ Age, data = wblake)
summary(ajuste)

Call:
lm(formula = Length ~ Age, data = wblake)

Residuals:
    Min      1Q  Median      3Q     Max 
-85.794 -19.499  -4.499  16.177  94.853 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  65.5272     3.1974   20.49   <2e-16 ***
Age          30.3239     0.6877   44.09   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 28.65 on 437 degrees of freedom
Multiple R-squared:  0.8165,    Adjusted R-squared:  0.8161 
F-statistic:  1944 on 1 and 437 DF,  p-value: < 2.2e-16

¿Resulta significativa la pendiente?
si
Resuma la bondad del ajuste con el R2:
el modelo explica un 82% de la variabilidad del largo

Superponga la recta estimada al gráfico de dispersión, y también las medias muestrales por grupos.

ggplot(wblake, aes(Age, Length))+
  geom_point()+
  geom_smooth(method = "lm", se = FALSE)+
  geom_boxplot(aes( color = as.factor(Age), fill=NULL, group = Age), alpha=0.5)+
  theme_minimal()+
  theme(legend.position = "none")

Halle el estimador de \(\sigma\) que proporciona el modelo lineal.

sigma(ajuste)

[1] 28.64585

¿A qué valor debiera parecerse?
se debe parecer al Residual standard error
¿Se parece?
si

sd(wblake$Length)

[1] 66.79134

¿Le parece que el ajuste obtenido por el modelo lineal es satisfactorio?
si

4. Obtenga intervalos de confianza de nivel 95% para la longitud media a edades 2, 4 y 6 años

new <- data.frame(Age = c(2,4,6))
#intervalo de confianza
predict.lm(ajuste,newdata = new,interval="confidence",level = 0.95)

       fit      lwr      upr
1 126.1749 122.1643 130.1856
2 186.8227 184.1217 189.5237
3 247.4705 243.8481 251.0929

¿Sería correcto obtener IC para la longitud media a los 9 años con este conjunto de datos?
implicaría realizar una extrapolación de los resultados.