Proy mod, mej

---
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(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)
```

Quitaremos observaciones influyentes previo a verificar los demás supuestos:

```{r}
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}
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. El gráfico Q-Q sugiere normalidad.

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.

Según la documentación, las actividades diarias (adla), habilidades motoras gruesas (grossmotor), y el índice de mobilidad (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 mobilidad del encuestado, por lo que sus valores se alinean bastante con los de mobilityind, dando problemas de multicolinealidad. Por ello, son removidas.

```{r}
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 significcativamente mayor a 1, por lo que 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=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)

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,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)))

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]],",",names[order[2]],", con",r2a[order[2]],",",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]]
```


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

# 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)
```

Comprobaremos los supuestos. Retirando primero observaciones influyentes:

```{r}
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. Asimismo, vemos que la varianza parece mantenerse constante. 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.

Utilizando mínimos cuadrados ponderados para la heterocedasticidad, y, posteriormente, retirando observaciones influyentes:

```{r}
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}
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}
estimation=estimation[-influcalc(smodelo4,estimation),]
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?