proyecto final 100% real no faik

---
title: "Proyecto Final Modelos I"
author: "Jorge Moreno, Aarón Calderón"
date: "2023-12-26"
geometry: "left=1.45cm,right=1.45cm,top=1.5cm,bottom=1.5cm"
output:
  pdf_document: default
  html_notebook: default
---

## Introducción

La fuerza máxima de agarre constituye un indicador esencial de la salud funcional, brindándonos valiosa información acerca de la capacidad física de un individuo. Este proyecto se centra en explorar y analizar dicha medida en países europeos, particularmente en España, utilizando la base de datos EasyShare, proveniente de la séptima encuesta de SHARE (Encuesta sobre Salud, Envejecimiento y Jubilación en Europa). Buscamos entender la relación entre la fuerza de agarre, factores demográficos e índices de limitación funcional. Este análisis proporcionará una comprensión integral de estas variables, permitiendo orientar estrategias de atención a la salud adaptadas a las necesidades específicas del envejecimiento y el bienestar de los españoles.

## Metodología

Considerando la base de datos EasyShare, lo primero que se realizó fue una limpieza de todas las observaciones, extrayendo únicamente las realizadas en España, y se eliminaron aquellas faltantes.
Mediante estadística descriptiva se analizó las variables categóricas y cuantitativas, y, a su vez, se examinó gráficamente la relación entre la medida de fuerza de agarre con cada una de las variables de interés.
Con el análisis realizado, se creó un modelo de regreso lineal con el fin de modelar la medida máxima de fuerza de agarre con base en las variables demográficas e índices de limitación funcional.
Del modelo completo, mediante la función StepAIC se escogió los mejores modelos para compararlos y elegir el mejor utilizando la media cuadrática del error de predicción. Con el mejor modelo para nuestro país se comprobó los supuestos y se proporcionó inferencias válidas que servirán de ayuda para los resultados. Por último, con todo lo obtenido, se concluye sobre la relación de las variables con respecto a la medida máxima de fuerza de agarre.

## Objetivo

Modelar la medida máxima de fuerza de agarre (maxgrip) en base de
variables demográficas (Demographics) y a índices de limitación funcional (Functional limitation
indices) para la población de España.

## Pregunta general de interés

¿Qué variables ayudan a explicar la medida máxima de fuerza de agarre?

## Proceso de modelación

### 1. Obtener el subconjunto de datos asignado a cada grupo y en base a las variables de estudio. Se debe realizar limpieza de la base de datos, si existe datos faltantes deben eliminar la fila de la observación. En la página 6 del pdf de la guía se encuentra la codificación que se usó para los datos faltantes.

No se incluye la variable religión, ya que esta fue solo preguntada en la wave 1. En el caso de España, no se realizó esa pregunta para ninguna observación. Asimismo, numeracy1 y numeracy2 no son incluidas porque tienen un gran número de valores faltantes.

```{r,echo=F}
#easyShare7_subset <- read.csv("C:/Users/USUARIO/Downloads/easyShare7_subset.csv")
easyShare7_subset <- read.csv("D:/Descargas/DATABASES/easyShare7_subset.csv")

share=subset(easyShare7_subset, country == 15,select = c(maxgrip,age,female,birth_country,citizenship,isced1997_r,mar_stat,age_partner,gender_partner,mobilityind,lgmuscle,adla,grossmotor,finemotor,iadlza,recall_1,recall_2))

find_aberrant_indices <- function(variable) {
  aberrant_values <- c(-3, -7, -9, -10, -11, -12, -13, -14, -15, -16)
  aberrant_indices <- which(variable %in% aberrant_values)
  return(aberrant_indices)
}
aberrant_indices_list <- lapply(share, find_aberrant_indices)
aberrant_indices <- unique(unlist(aberrant_indices_list))
share_cleaned <- share[-aberrant_indices, ]
```

### 2. Realizar una estadística descriptiva de la muestra con base en las variables de estudio.

```{r,echo=F}
library(knitr)

a=as.data.frame(round(rbind(summary(share_cleaned$maxgrip),summary(share_cleaned$age),summary(share_cleaned$age_partner),summary(share_cleaned$recall_1),summary(share_cleaned$recall_2)),2))
rownames(a)=c("Fuerza de agarre máxima","Edad","Edad de la pareja","Recordar palabras 1","Recordar palabras 2")
colnames(a)=c("Mínimo","Cuartil 1","Mediana","Media","Cuartil 3","Máximo")

kable(a,caption="Medidas de tendencia central")
```

En la fuerza de agarre máxima, considerando que el máximo es 100, vemos que son valores bastante bajos, siendo el tercer cuartil 37. Esto está explicado por la siguiente variable, edad, donde vemos que la mayoría de observaciones son de personas de la tercera edad. La edad de la pareja es bastante similar a la variable edad.

Finalmente, las dos variables de recordar palabras dan el número de palabras que el encuestado puede recordar por primera vez (Recordar 1) o después de un tiempo (Recordar 2). En general, recuerdan 5 palabras por primera vez, y 3 después de un tiempo.

```{r,echo=F}
genero=t(as.data.frame(table(share_cleaned$female))[-1])
por1=round(genero/sum(genero)*100,2)
por1=sprintf("%.2f%%",por1)
genero=rbind(genero,por1)
colnames(genero)=c("Masculino","Femenino")
rownames(genero)=c("Género","Porcentaje")
kable(genero,caption="Género del encuestado")
```

```{r,echo=F}
pais=table(share_cleaned$birth_country)
noespana=sum(pais[-27])
espana=pais[27]
paiss=data.frame(noespana,espana)

por2=round(paiss/sum(paiss)*100,2)
por2=sprintf("%.2f%%",por2)
paiss=rbind(paiss,por2)
colnames(paiss)=c("Fuera de España","España")
rownames(paiss)=c("País natal","Porcentaje")


kable(paiss,caption="País de origen del encuestado.")
```

La gran mayoría de los encuestados son españoles, siendo apenas un 5% extranjeros.

```{r,echo=F}
ciu=table(share_cleaned$citizenship)
noesp=sum(ciu[-18])
esp=ciu[18]
ciuu=data.frame(esp,noesp)
por3=round(ciuu/sum(ciuu)*100,2)
por3=sprintf("%.2f%%",por3)
ciuu=rbind(ciuu,por3)
colnames(ciuu)=c("Sí","No")
rownames(ciuu)=c("Posesión","Porcentaje")
kable(ciuu,caption="Posesión de la ciudadanía española")
```

Vemos que casi todos poseen la nacionalidad española.

```{r,echo=F}
ed=t(as.data.frame(table(share_cleaned$isced1997_r))[-1])
por4=round(ed/sum(ed)*100,2)
por4=sprintf("%.2f%%",por4)
ed=rbind(ed,por4)
colnames(ed)=c("Ninguno","1","2","3","4","5","6","En curso","Otro")
rownames(ed)=c("Nivel","Porcentaje")

kable(ed,caption="Nivel educativo del encuestado")
```

Vemos que la mayoría de observaciones no tuvieron estudios, tuvieron estudios primarios, o inicios de los secundarios. Parece haber un salto de personas que llegaron a la primera parte de estudios terciarios, en contraste con los que no terminaron los secundarios.


```{r,echo=F}
casao=t(as.data.frame(table(share_cleaned$mar_stat))[-1])
por5=round(casao/sum(casao)*100,2)
por5=sprintf("%.2f%%",por5)
casao=rbind(casao,por5)
casao=cbind(c(0,0),casao)
colnames(casao)=c("Ninguno","Casado, mismo hogar","Unión registrada","Casado, no mismo hogar","Nunca casado","Divorciado","Viudo")
rownames(casao)=c("Estado Civil","Porcentaje")

kable(casao,caption="Estado civil del encuestado")
```

La gran mayoría de observaciones están casadas y viven en el mismo hogar.

```{r,echo=F}
genero2=t(as.data.frame(table(share_cleaned$gender_partner))[-1])
por6=sprintf("%.2f%%",round(genero2/sum(genero2)*100,2))
genero2=rbind(genero2,por6)
colnames(genero2)=c("Masculino","Femenino")
rownames(genero2)=c("Género","Porcentaje")
kable(genero2,caption="Género de la pareja del encuestado")
```

Valores muy similares a los del género del encuestado.
A continuación, se presentará la información de los índices de limitación funcional. Recuérdese que a mayor índice, mayor dificultad en los parámetros evaluados.

```{r,echo=F}
mob2=t(as.data.frame(table(share_cleaned$mobilityind))[-1])
por7=sprintf("%.2f%%",round(mob2/sum(mob2)*100,2))
mobi2=rbind(mob2,por7)
colnames(mobi2)=c("0","1","2","3","4")
rownames(mobi2)=c("Índice","Porcentaje")
kable(mobi2,caption="Índice de movilidad del encuestado")

lgm2=t(as.data.frame(table(share_cleaned$lgmuscle))[-1])
por8=sprintf("%.2f%%",round(lgm2/sum(lgm2)*100,2))
lgmu2=rbind(lgm2,por8)
colnames(lgmu2)=c("0","1","2","3","4")
rownames(lgmu2)=c("Índice","Porcentaje")
kable(lgmu2,caption="Índice del músculo grande del encuestado")

adl2=t(as.data.frame(table(share_cleaned$adla))[-1])
por9=sprintf("%.2f%%",round(adl2/sum(adl2)*100,2))
adla2=rbind(adl2,por9)
colnames(adla2)=c("0","1","2","3","4","5")
rownames(adla2)=c("Índice","Porcentaje")
kable(adla2,caption="Índice de actividades diarias del encuestado")

gro2=t(as.data.frame(table(share_cleaned$grossmotor))[-1])
por10=sprintf("%.2f%%",round(gro2/sum(gro2)*100,2))
gros2=rbind(gro2,por10)
colnames(gros2)=c("0","1","2","3","4")
rownames(gros2)=c("Índice","Porcentaje")
kable(gros2,caption="Índice de las habilidades motoras gruesas del encuestado")

fin2=t(as.data.frame(table(share_cleaned$finemotor))[-1])
por11=sprintf("%.2f%%",round(fin2/sum(fin2)*100,2))
fine2=rbind(fin2,por11)
colnames(fine2)=c("0","1","2","3")
rownames(fine2)=c("Índice","Porcentaje")
kable(fine2,caption="Índice de las habilidades motoras finas del encuestado")

iad2=t(as.data.frame(table(share_cleaned$iadlza))[-1])
por12=sprintf("%.2f%%",round(iad2/sum(iad2)*100,2))
iadl2=rbind(iad2,por12)
colnames(iadl2)=c("0","1","2","3","4","5")
rownames(iadl2)=c("Índice","Porcentaje")
kable(iadl2,caption="Índice de las actividades instrumentales diarias del encuestado")
```

Vemos que los encuestados, en su mayoría, no presentan anomalías con las limitaciones funcionales, siendo la mayor excepción el índice de músculo grande (actividades que requieren músculos de las piernas), donde hay una cantidad importante de individuos que presentan una cierta dificultad con ello.

### 3. Realizar una exploración gráfica de la relación entre la variable de respuesta y las variables explicativas.

```{r,echo=F}
ygraf=share_cleaned$maxgrip
par(mfrow=c(3,4),mai=c(0.50,0.50,0.50,0.50))
plot(ygraf,x=share_cleaned$age,ylab="Agarre",xlab="Edad")
plot(ygraf,x=share_cleaned$female,ylab="Agarre",xlab="Género")
plot(ygraf,x=share_cleaned$birth_country,ylab="Agarre",xlab="País nac.")
#plot(ygraf,x=share_cleaned$citizenship,ylab="Agarre",xlab="Ciudadanía")
#plot(ygraf,x=share_cleaned$isced1997_r,ylab="Agarre",xlab="Educación")
#plot(ygraf,x=share_cleaned$isced1997_r,xlim=c(-1,7),ylab="Agarre",xlab="Educación (limitada)")
#plot(ygraf,x=share_cleaned$mar_stat,ylab="Agarre",xlab="Estado civil")
plot(ygraf,x=share_cleaned$age_partner,ylab="Agarre",xlab="Edad pareja")
plot(ygraf,x=share_cleaned$gender_partner,ylab="Agarre",xlab="Género pareja")
plot(ygraf,x=share_cleaned$mobilityind,ylab="Agarre",xlab="Índ. movilidad")
plot(ygraf,x=share_cleaned$lgmuscle,ylab="Agarre",xlab="Índ. músculo grande")
plot(ygraf,x=share_cleaned$adla,ylab="Agarre",xlab="Índ. actividades diarias")
plot(ygraf,x=share_cleaned$grossmotor,ylab="Agarre",xlab="Índ. hab. motoras gruesas")
plot(ygraf,x=share_cleaned$finemotor,ylab="Agarre",xlab="Índ. hab. motoras finas")
plot(ygraf,x=share_cleaned$iadlza,ylab="Agarre",xlab="Índ. actividades instrumentales")
plot(ygraf,x=share_cleaned$recall_1,ylab="Agarre",xlab="Recordar palabras 1")
plot(ygraf,x=share_cleaned$recall_2,ylab="Agarre",xlab="Recordar palabras 2")
```



### 4. Proponer un modelo completo en base a la exploración gráfica. Comprobar los supuestos del modelo y en caso de ser necesario, corregir. Realizar selección de variables. Proponer 3 mejores modelos.

Se explicará las variables que se considerará, con base en la exploración gráfica.
Observemos que a medida que avanza la edad, el agarre disminuye.  En el caso de género, el género masculino tiene mayor fuerza de agarre. Sucede una situación similar con la edad de la pareja, y género de la pareja, donde los individuos con parejas femeninas tienen mayor agarre.

La movilidad, el músculo, actividades diarias, habilidades motoras gruesas y finas, y actividades dirarias instrumentales guardan la misma relación: a medida que aumentan (esto es, en la vida real, empeoran), el agarre disminuye. Esto es mucho más sutil en músculo.

El agarre parece acrecentar cuando crece la habilidad de recordar palabras 1. En cuanto al resto de variables: país de nacimiento, ciudadanía, educación, estado civil, y habilidad de recordar palabras 2 no parecen guardar relación con el agarre, ya sea por la gráfica, o porque no tienen relación aparente. Por tanto, estas no se consideran para el modelo. La estructura del modelo sería:

$$
\text{fuerza}=\beta_{0}+\beta_{1}\text{edad}+\beta_{2}\text{sexo}+\beta_{3}\text{movilidad1}+\beta_{4}\text{movilidad2}+\beta_{5}\text{movilidad3}+\beta_{6}\text{movilidad4}+\beta_{7}\text{palabras}_{1}
$$

Como se tiene 2977 observaciones, se usará 992 observaciones para los conjuntos de entrenamiento, validación y estimación.

```{r,echo=F,message=F,warning=F}
suppressWarnings(library(dplyr,quietly=T))
copy=share_cleaned[,c(-4,-5,-6,-7,-17)]

copy <- copy %>%
  mutate(female = ifelse(female == 0, "Male", "Female"))

copy <- copy %>%
  mutate(gender_partner = ifelse(gender_partner == 0, "Male", "Female"))

copy$female=factor(copy$female)
copy$gender_partner=factor(copy$gender_partner)
copy$mobilityind=factor(copy$mobilityind)
copy$lgmuscle=factor(copy$lgmuscle)
copy$adla=factor(copy$adla)
copy$grossmotor=factor(copy$grossmotor)
copy$finemotor=factor(copy$finemotor)
copy$iadlza=factor(copy$iadlza)

n=nrow(copy)

set.seed(1291)
samp=sample(1:2977, size= 992  ,replace = F)
set.seed(1171)
samp2=sample(c(1:2977)[-samp], size= 992  ,replace = F)
training=copy[samp,]
validation=copy[-c(samp,samp2),]
estimation=copy[samp2,]


modelo_com=lm(maxgrip~age+female+age_partner+
                gender_partner+mobilityind+lgmuscle+adla+grossmotor+finemotor+
                iadlza+recall_1,data = training)
```

Después de quitar observaciones influyentes previo a verificar los demás supuestos:

```{r,echo=F}
calcular_anomalias <- function(modelo, dataset) {
  dffits <- dffits(modelo)
  cook_dist <- cooks.distance(modelo)
  dfbetas_mat <- dfbetas(modelo)

  c_valor <- 2*sqrt(length(coef(modelo)) / nrow(dataset))
  d_valor <- 2 / sqrt(nrow(dataset))

  indices_dffits <- which(abs(dffits) > c_valor)
  indices_cook <- which(cook_dist > 1)
  indices_influyentes <- numeric(0)
  for (i in 1:ncol(dfbetas_mat)) {
    indices_influyentes <- union(indices_influyentes, which(abs(dfbetas_mat[,i]) > d_valor))
  }
  indices_anomalos <- unique(c(indices_dffits, indices_cook, indices_influyentes))
  return(indices_anomalos)
}

modelo_com=lm(maxgrip~age+female+age_partner+
                gender_partner+mobilityind+lgmuscle+adla+grossmotor+finemotor+
                iadlza+recall_1,data = training)

training=training[-calcular_anomalias(modelo_com,training),]

modelo_com=lm(maxgrip~age+female+age_partner+
                gender_partner+mobilityind+lgmuscle+adla+grossmotor+finemotor+
                iadlza+recall_1,data = training)
```


```{r,echo=F,warning=F,message=F,out.width="85%",fig.align='center'}
par(mfrow=c(2,2))
plot(modelo_com,cex=0.5)
```

```{r,echo=F,out.width="50%",fig.align='center'}
plot(rstandard(modelo_com))
```


Podemos ver que se cumple el supuesto de independencia, ya que no notamos la existencia de algún patrón, hay normalidad, no cumple el supuesto de homocedasticidad ya que observamos que la variabilidad de los datos aumenta, se puede apreciar linealidad y hay presencia de valores aberrantes. El gráfico Q-Q sugiere normalidad.

Se utilizará mínimos cuadrados ponderados para corregir la hetercedasticidad. Los residuos ponderados después de este cambio:

```{r,echo=F,out.width="50%",fig.align='center'}
mod=lm(abs(modelo_com$residuals)~modelo_com$fitted.values)
w=1/(mod$fitted.values^2)
formula1=formula(modelo_com)
modelo_1=lm(formula1,data = training,weights =w)
respond=sqrt(w)*modelo_1$residuals

plot(modelo_1$fitted.values,respond)
```

Corregiremos multicolinealidad. Para ello, se quitará variables que tengan mucha correlación entre ellas.

Es lógico que Age y Age_partner, al igual que gender y gender_partner estén muy correlacionadas, pues las parejas son de edades similares, y, frecuentemente, son del género opuesto del encuestado. Por ello, se quitará Age_partner y gender_partner.

Según la documentación, las actividades diarias (adla), habilidades motoras gruesas (grossmotor), y el índice de movilidad (mobilityind) miden habilidad para caminar, bañarse o subir escaleras, por lo que sus valores están correlacionados, y producirán multicolinealidad. Tambien, el músculo grande (lgmuscle), está relacionado con actividades que involucran al estado de las piernas. Por tal motivo, se conservará solo mobilityind.

Asimismo, las habilidades motoras finas (finemotor) y las actividades instrumentales diarias (iadlza), tienen actividades similares, como recoger objetos pequeños, preparar o cortar comida. De igual forma, ambas pretenden medir la movilidad del encuestado, por lo que sus valores se alinean bastante con los de mobilityind, dando problemas de multicolinealidad. Por ello, son removidas. Calculando los VIF después de estos cambios, junto con su media:

```{r,message=F,warning=F,echo=F}
suppressWarnings(library(car))
modelo0=update(modelo_1,~.-age_partner-finemotor-grossmotor-gender_partner-adla-lgmuscle-iadlza)
vif(modelo0)
mean(vif(modelo0)[,3])
```

Se considerará que la media no es significativamente mayor a 1, por lo que no hay multicolinealidad.

Para escoger los mejores modelos, se utilizará stepAIC, donde se incluirá modelos con efectos cuadráticos y efectos de interacción, usando los criterios de AIC y BIC. Luego, se los seleccionará a partir de su $R^{2}$ ajustado.

```{r,echo=F,message=F,warning=F}
library(MASS)
quad=function(vec){
  return((vec-mean(vec))^2)
}
copy2=training[,c(1,2,3,6,12)]
copy2$age2=quad(training$age)
copy2$recall_12=quad(training$recall_1)

form=formula(modelo0)
modelo01=lm(form,data=copy2)
uppermod1=lm(maxgrip~.,data=copy2)

smodelo1=stepAIC(modelo01,scope=list(upper=uppermod1),direction="both",trace=F)
smodelo2=stepAIC(modelo01,scope=list(upper=uppermod1),direction="both",trace=F,k=log(nrow(training)))

copy3=training[,c(1,2,3,6,12)]
uppermod2=lm(maxgrip~.^2,data=copy3)
smodelo3=stepAIC(modelo0,scope=list(upper=uppermod2),direction="both",trace=F)
smodelo4=stepAIC(modelo0,scope=list(upper=uppermod2),direction="both",trace=F,k=log(nrow(training)))

topmod=function(model_list,names) {
  r2a=sapply(model_list, function(model) summary(model)$adj.r.squared)
  order=sort(r2a,decreasing=T,index.return=T)$ix
  cat("Modelos con mayor R2 ajustado:",names[order[1]],", con",r2a[order[1]],"\n")
  cat(names[order[2]],", con",r2a[order[2]],"\n")
  cat(names[order[3]],", con",r2a[order[3]])
  return(c(model_list[order[1]],model_list[order[2]],model_list[order[3]]))
}

mmods=topmod(list(modelo0,smodelo1,smodelo2,smodelo3,smodelo4),c("Modelo original","Cuadrático AIC","Cuadrático BIC","Interacción AIC","Interacción BIC"))
mmod1=mmods[1][[1]]
mmod2=mmods[2][[1]]
mmod3=mmods[3][[1]]
```
Los 3 mejores modelos en el mismo orden mencionado son:

$$
\begin{aligned}
\text{maxgrip}&=39.48-0.26age+26.07female+11.87mobilityind1-6.62mobilityind2+5.06mobilityind3+90.33mobilityind4+\\
&0.38recall_1-0.16agefemale-0.2agemobilityind1+0.05agemobilityind2-0.13agemobilityind3-1.18agemobilityind4\\
\text{maxgrip}&=41.74-0.3age+25.26female-2.86mobilityind1-3.04mobilityind2-5.03mobilityind3-6.15mobilityind4\\
&+0.40recall_1-0.14agefemale\\
\text{maxgrip}&=44.42-0.34age+15.14female-2.76mobilityind1-3.13mobilityind2-4.91mobilityind3-7.24mobilityind4\\
&+0.42recall_1
\end{aligned}
$$

### 5. Validar los modelos y escoger uno. Comprobar supuestos

Calculamos la media cuadrática del error de predicción MCEP de los 3 mejores modelos y los comparamos con la media cuadrática del error MCE.  En el mismo orden para los modelos listados anteriormente:

```{r,echo=F}
validation$age2=quad(validation$age)
validation$recall_12=quad(validation$recall_1)

# Mejor modelo 1
msepAIC_val=mean((validation$maxgrip-predict(mmod1,validation))^2)
mce1=sigma(mmod1)^2
(dif1=abs(msepAIC_val-mce1))

#  Mejor modelo 2
msepBIC_val=mean((validation$maxgrip-predict(mmod2,validation))^2)
mce2=sigma(mmod2)^2
(dif2=abs(msepBIC_val-mce2))

# Mejor modelo 3
msep_val=mean((validation$maxgrip-predict(mmod3,validation))^2)
mce3=sigma(mmod3)^2
(dif3=abs(msep_val-mce3))

modelovalid=mmods[which.min(c(dif1,dif2,dif3))][[1]]
vmodelo=lm(formula(modelovalid),data=validation)
```

El mejor modelo es el del BIC. omprobaremos los supuestos. Retirando primero observaciones influyentes:

```{r,echo=F,warning=F,message=F,out.width="85%",fig.align='center'}
validation=validation[-calcular_anomalias(vmodelo,validation),]

vmodelo=lm(formula(modelovalid),data=validation)
par(mfrow=c(2,2))
plot(vmodelo)
```

Ya que no hay ningún patrón en la gráfica de valores ajustados y residuos, hay linealidad. La variabilidad de los datos no es constante, por lo que puede haber problemas de heterocedasticidad. El gráfico Q-Q nos dice que sí hay normalidad.

Se usará mínimos cuadrados ponderados con el fin de resolver la heterocedasticidad:

```{r,echo=F,out.width="50%",fig.align='center'}
tmod=lm(abs(vmodelo$residuals)~vmodelo$fitted.values)
ws=1/(tmod$fitted.values^2)
vmodelo1=lm(formula(modelovalid),data = validation,weights =ws)
pondres=sqrt(ws)*vmodelo1$residuals

plot(vmodelo1$fitted.values,pondres)
```

Verificando multicolinealidad:

```{r,echo=F}
vif(vmodelo,type="predictor")
mean(vif(vmodelo,type="predictor")[,3])
```

### 6. Estimar el modelo. Proporcionar *inferencias válidas* por medio de pruebas de significancia.

```{r,echo=F,warning=F,message=F,out.width="85%",fig.align='center'}
modest=lm(formula(vmodelo1),data=estimation)
estimation=estimation[-calcular_anomalias(modest,estimation),]
modest=lm(formula(vmodelo1),data=estimation)
par(mfrow=c(2,2))
plot(modest)
```

Corrigiendo heterocedasticidad:

```{r,echo=F,out.width="50%",fig.align='center'}
emod=lm(abs(modest$residuals)~modest$fitted.values)
ws2=1/(emod$fitted.values^2)
modest1=lm(formula(modest),data = estimation,weights=ws2)
pondres2=sqrt(ws2)*modest1$residuals

plot(modest1$fitted.values,pondres2)
summary(modest1)
```

Si usamos un nivel de significancia del 5%, las variables edad, género, índice de movilidad, la habilidad de recordar palabras y el efecto interactivo entre edad y género tienen un valor p menor que dicho nivel. Por lo tanto, se rechaza la hipótesis nula que indica que los coeficientes de estas variables son 0, y se concluye que son significativas.

$\beta_{0}=36.95$ es el agarre promedio para una mujer recién nacida, con índice de movilidad 0, y que no recuerda ninguna palabra de las que se pidió recordar.

$\beta_{1}=-0.23$ es el efecto adicional en el agarre promedio por cada año de edad adicional para una mujer con cualquier índice de movilidad, suponiendo que el número de palabras recordadas se mantiene constante.

$\beta_{2}=31.34$ es el efecto adicional en el agarre promedio para un hombre recién nacido, con índice de movilidad 0, y que no recuerda ninguna palabra de las que se le pidió recordar.

$\beta_{3}=-1.48,\beta_{4}=-3.79\text{ y }\beta_{5}=-6.41$ son la diferencia en el agarre promedio de una persona con índice de movilidad 0 y otra con índice de movilidad 1, 2 o 3, respectivamente, del mismo sexo, y suponiendo que tienen la misma edad y recuerdan el mismo número de palabras.

Para nuestro nivel de significancia del 5%, no existe evidencia para rechazar la hipótesis nula que indica que $\beta_{6}=0$, según la prueba t. Por lo tanto, no se interpreta.

$\beta_{7}=0.6$ es el efecto adicional en el agarre promedio por cada palabra adicional recordada para un hombre o mujer cuando el índice de movilidad es el mismo y la edad se mantiene constante.

$\beta_{8}=-0.24$ por cada año adicional, y cuando el índice de movilidad y palabras recordadas es constante, el agarre promedio cambia en $-0.24$ kg. más que lo cambia en las mujeres. Así, por cada año adicional, con movilidad y palabras constantes, el cambio en el agarre promedio de un hombre es de $\beta_{8}+\beta_{1}=-0.47$.

A continuación, se obtendrá intervalos de confianza para los coeficientes anteriormente estimados. Esto es, con un 95% de confianza, el verdadero valor del coeficiente estará dentro del siguiente rango:

```{r,echo=F}
confint(modest1,level = 0.95)
```

### 7. Concluir con base en la pregunta de investigación: ¿qué variables muestran una relación positiva con la respuesta? ¿Qué variables muestran una relación negativa con la respuesta?

Es bueno recalcar, la importancia de definir el país en el que uno quiere modelar para extraer y reducir la base de datos, de igual manera, eliminar todas aquellas observaciones aberrantes e influyentes que cuente la base de datos, de modo que se puedan obtener conclusiones más precisas.

Para concluir, en caso de España no se tomó en las encuestas realizadas la variable religión y presentaron muy pocas observaciones las variables  numeracy1 y numeracy2. Una vez analizadas todas las variables de interés, se eliminó algunas variables por su relación con otras y su poca importancia dentro del objetivo del estudio.

Luego de toda la metodología aplicada, las variables explicativas que mejor explican la medida máxima de fuerza de agarre son la edad, sexo, índice de movilidad, la habilidad de recordar palabras y el efecto interactivo entre edad y género. La edad presenta una relación negativa con la respuesta, esto es, a mayor edad, menor fuerza de agarre en promedio.

Por el contrario, la capacidad de recordar palabras tiene una relación positiva: mientras más palabras se recuerde, mayor fuerza de agarre se tiene en promedio.

Para género y el índice de movilidad, el ser hombre tiene un efecto positivo en la fuerza de agarre, esto es, la incrementa, en comparación con el ser mujer. Por otro lado, el tener un índice de movilidad distinto de 0 tiene un efecto negativo en la variable de respuesta, en comparación con un índice de movilidad de 0.