proy mod corr

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

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 mobilidad 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}
par(mfrow=c(3,3),mai=c(0.55,0.55,0.55,0.55))
ygraf=share_cleaned$maxgrip
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. mobilidad")
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 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"))

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)
n/3
set.seed(2024)
samp=sample(1:2977, size= 992  ,replace = F)
set.seed(2025)
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)
```

```{r,echo=F}
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,echo=F}
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.

Ahora, las habilidades motoras y las actividades diarias están correlacionadas, pues se usan las primeras para desarrollar estas últimas. Por ello, se eliminará adla.

Finalmente, las habilidads motoras y el índice de mobilidad, al igual que con el músculo grande, también están correlacionadas, y usarlas podría darnos información repetida. Por ello, se usará solo la habilidad motora gruesa, y el índice de mobilidad. Se quitará grossmotor, y lgmuscle.


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

No hay multicolinealidad. Se procederá a quitar observaciones influyentes:

```{r}
influcalc <- function(modelo, dataset) {

  dffits <- dffits(modelo)
  cook_dist <- cooks.distance(modelo)

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

  indices_dffits <- which(abs(dffits) > c_valor)
  indices_cook <- which(cook_dist > 1)

  indices_anomalos <- unique(c(indices_dffits, indices_cook))
  return(indices_anomalos)
}

training=training[-influcalc(modelo0,training),]

noinflu=formula(modelo0)
modelo0=lm(noinflu,data=training)
```

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=training[,c(1,2,3,6,9,11,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)

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=training[,c(1,2,3,6,9,11,12)]
uppermod2=lm(maxgrip~.^2,data=copy2)
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 interacciones por BIC, y el modelo con efectos cuadráticos por AIC.

### 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}
validation$age2=quad(validation$age)
validation$recall_12=quad(validation$recall_1)

# 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(smodelo1)^2
(dif2=abs(msepBIC_val-mce2))

# modelo con efectos cuadráticos por AIC
msep_val=mean((validation$maxgrip-predict(smodelo1,validation))^2)
mce3=sigma(modelo0)^2
(dif3=abs(msep_val-mce3))
```
El modelo con efectos de interacción escogido por el BIC 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(smodelo4)
```

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}
modelo3_sinpeso=lm(maxgrip ~ age + female + mobilityind +
    recall_1 + age:female, data = training)
vif(modelo3_sinpeso,type="predictor")
```

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