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

## Introducción

SHARE, la Encuesta sobre Salud, Envejecimiento y Jubilación en Europa, es una infraestructura de
investigación para estudiar los efectos de las políticas sanitarias, sociales, económicas y medioambientales a lo
largo de la vida de los ciudadanos europeos y más allá.

Para este proyecto utilizaremos EasyShare que es una base de datos para uso académico que condensa parte
de la información de SHARE de 7 encuestas (waves).

## 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.

```{r}
#easyShare7_subset <- read.csv("C:/Users/USUARIO/Downloads/easyShare7_subset.csv")
easyShare7_subset <- read.csv("D:/Descargas/easyShare7_subset.csv")
attach(easyShare7_subset)
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))
# No fue incluída q34_re (religion) ya que esta solo fue incluida en la wave 1
# En caso de España, no se realizó la pregunta de religión para ninguna observación.
# No se incluyen orienti, numeracy_1 y numeracy_2 por tener un gran número de valores aberrantes.
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)
detach(easyShare7_subset)
attach(share_cleaned)

s1=summary(maxgrip)
s2=summary(age)
s3=summary(age_partner)
s4=summary(mobilityind)
s5=summary(lgmuscle)
s6=summary(adla)
s7=summary(grossmotor)
s8=summary(finemotor)
s9=summary(iadlza)
s10=summary(recall_1)
s11=summary(recall_2)
a=as.data.frame(round(rbind(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11),2))
rownames(a)=c("Fuerza de agarre máxima","Edad","Edad de la pareja","Índice de mobilidad","Índice muscular","Índice de actividades diarias","Habilidades motoras gruesas","Habilidades motoras finas","Índice de actividades diarias instrumentales","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. Para las demás variables, con excepción de las de recordar palabras, vemos valores extremadamente bajos. Recordemos que estas toman valores de 0 hasta 3, 4 o 5, dependiendo de la variable, y mientras más alto sea el valor, mayor es la dificultad que tienen en la medida examinada.

Observemos que, en general, no hay muchos problemas en estas áreas, con la excepción del índice muscular, que es un poco más elevado.

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(female))[-1])
por1=round(genero/sum(genero)*100,2)
genero=rbind(genero,por1)
colnames(genero)=c("Masculino","Femenino")
rownames(genero)=c("Género del encuestado","Porcentaje")
kable(genero)
```

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

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


kable(paiss)
```

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

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

Vemos que casi todos poseen la nacionalidad española.


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

kable(ed)
```

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(mar_stat))[-1])
por5=round(casao/sum(casao)*100,2)
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)
```

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

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

Valores muy similares a los del género del encuestado.

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

```{r,echo=F}
par(mfrow=c(3,3),mai=c(0.55,0.55,0.55,0.55))
plot(y=maxgrip,x=age)
plot(y=maxgrip,x=female)
plot(y=maxgrip,x=birth_country)
plot(y=maxgrip,x=citizenship)

plot(y=maxgrip,x=isced1997_r)
plot(y=maxgrip,x=isced1997_r,xlim=c(-1,7))
plot(y=maxgrip,x=mar_stat)
plot(y=maxgrip,x=age_partner)
plot(y=maxgrip,x=gender_partner)

plot(y=maxgrip,x=mobilityind)
plot(y=maxgrip,x=lgmuscle)
plot(y=maxgrip,x=adla)
plot(y=maxgrip,x=grossmotor)

plot(y=maxgrip,x=finemotor)
plot(y=maxgrip,x=iadlza)
plot(y=maxgrip,x=recall_1)
plot(y=maxgrip,x=recall_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 mobilidad, 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 aagarre. Por tanto, estas no se consideran para el modelo.

```{r,echo=F}
suppressWarnings(library(dplyr))
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"))
```


```{r}
detach(share_cleaned)
suppressWarnings(attach(copy))
n=nrow(copy)
n/3
set.seed(451)
samp=sample(1:2977, size= 992  ,replace = F)
set.seed(124)
samp2=sample(c(1:2977)[-samp], size= 992  ,replace = F)
training=copy[samp,]
validation=copy[-c(samp,samp2),]
estimation=copy[samp2,]

female=factor(female)
gender_partner=factor(gender_partner)


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

```{r}
par(mfrow=c(2,2))
plot(modelo_com)

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. Se utilizará mínimos cuadrados ponderados para corregir la hetercedasticidad.

```{r}
mod=lm(abs(modelo_com$residuals)~modelo_com$fitted.values)
w=1/(mod$fitted.values^2)
#birth_country+citizenship+mar_stat
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:

```{r}
suppressWarnings(library(car))
vif(modelo_1)
cor(copy[,c(1,2,4,6:12)])
```

Age y Age_partner, adla y finemotor, grossmotor y adla, lgmuscle y mobilityidnex y grossmotor y mobilityindex guardan fuertes correlaciones entre sí. Los géneros de las parejas también tienen una gran correlación. Se eliminará: age_partner, grossmotor, adla, lgmuscle y gender_partner:

```{r}
modelo0=update(modelo_1,~.-age_partner-grossmotor-gender_partner-adla-lgmuscle)
vif(modelo0)
mean(vif(modelo0))
```

Ya no hay multicolinealidad. Para escoger los mejores modelos, se utilizará stepAIC, donde el modelo máximo será: el modelo completo, el modelo con efectos cuadráticos, y el modelo con efectos de interacción.

```{r}
quad=function(vec){
  return((vec-mean(vec))^2)
}

copy2=copy[,c(1,2,3,6,10,11,12)]
copy2$age2=quad(age)
copy2$mobilityind2=quad(mobilityind)
copy2$finemotor2=quad(finemotor)
copy2$iadlza2=quad(iadlza)

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

library(MASS)
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)))

copy2=copy[,c(1,2,3,6,10,11,12)]
uppermod2=lm(maxgrip~.^2,data=copy)
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)))
```


Los modelos con un mayor $R^{2}$ ajustado, y que pasarán a la siguiente fase son: el modelo con efectos de interacción por AIC, el modelo con efectos de interacción por BIC, y el modelo original.

### 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.


```{r}
# modelo con efectos de interacción por AIC
msepAIC_val=mean((validation$maxgrip-predict(smodelo3,validation))^2)
mce1=sigma(smodelo3)^2
(dif1=abs(msepAIC_val-mce1))

#  modelo con efectos de interacción por BIC
msepBIC_val=mean((validation$maxgrip-predict(smodelo4,validation))^2)
mce2=sigma(smodelo4)^2
(dif2=abs(msepBIC_val-mce2))

# modelo original
msep_val=mean((validation$maxgrip-predict(modelo0,validation))^2)
mce3=sigma(modelo0)^2
(dif3=abs(msep_val-mce3))
```
El modelo con efectos de interacción escogido por el AIC es el que presenta menor discrepancia entre el MSE y el MSEP, por tanto, ese será el modelo escogido. Comprobaremos los supuestos:

```{r}
par(mfrow=c(2,2))
plot(smodelo3)
```

Ya que no hay ningún patrón en la gráfica de valores ajustados y residuos, hay linealidad. Asimismo, vemos que la varianza parece mantenerse constante, por lo que suponemos que hay homocedasticidad.
El criterio de normalidad parece cumplirse.

```{r}
vif(smodelo3,type="predictor")
```

Averiguar.

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

```{r}
formula999=formula(smodelo3)
modelofinal=lm(formula999,data=estimation)
summary(modelofinal)
anova(modelofinal)
```

Si usamos un nivel de significancia del 5%, las variables edad, género, índice de mobilidad, actividades diarias instrumentales, 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.

```{r}
#VERIFICAR CÓMO DEBE HACERSE PARA MÍNIMOS CUADRADOS PONDERADOS!!
confint(modelofinal)
```


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