Taller 6 modelos

---
title: "Taller 6_Modelos I"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# Set so that long lines in R will be wrapped:
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=70),tidy=TRUE)
```

# Unidad 8: Construcción del modelo de regresión lineal simple.


## Dataset: SMA

Considere los datos sobre 141 áreas metropolitanas estándar en los Estados Unidos (sma.txt). Los datos proporcionan información sobre 11 variables para cada área para el período 1976-1977. Las zonas se han dividido en 4 regiones geográficas: 1=Noreste, 2=Norte-Centro, 3=Sur,
4=Oeste. Las variables proporcionadas son las siguientes:


Área: tamaño de la tierra en millas cuadradas
Población total: estimada en miles
Porcentaje ciudad:  porcentaje de la población en el centro de la ciudad / ciudades
Porcentaje senior: porcentaje de personas mayores de $\leq$ 65 años
Médicos: número de médicos profesionalmente activos
Camas de hospital: número total de camas de hospital
Graduados:  porcentaje de adultos que terminaron la escuela secundaria
Laboral: número de personas en la fuerza laboral en miles
Ingresos: ingresos totales en 1976 en millones de dólares
Delitos: número de delitos graves
Región: región geográfica según el Censo de los Estados Unidos

La variable respuesta que estamos interesados en modelar es la tasa de delitos por cada mil habitantes (Delitos/Población total).

```{r}
sma <- read.table("D:/Descargas/DATABASES/sma.txt", quote="\"", comment.char="")
sma=sma[,-1]
names(sma) <- c('Area', 'PobT', 'PCiudad', 'PSenior', 'Medicos', 'CHosp', 'Grad', 'Laboral', 'Ingresos', 'Delitos', 'Region')
head(sma)
#summary(sma)


# Variable de respuesta
sma$TDelitos <- sma$Delitos/sma$PobT
summary(sma$TDelitos)
hist(sma$TDelitos, xlab = "Tasa de delitos", main = "Histograma de la tasa de delitos")
boxplot(sma$TDelitos)
```

## 1. Divida aleatoriamente el set de datos en tres partes, una para la construcción del modelo, otra para la validación y finalmente otra para la estimación final. Compruebe que el tamaño del set de construcción del modelo y la estimación final cumpla con el requerimiento mínimo (n $>$ 5P, recordando que para la variable región, se estimaría 3 parámetros). Utilice como P el número de coeficientes en un modelo aditivo de primer orden incluyendo todas las variables explicativas.


```{r}
sma$Region_n=sma$Region
sma$Region_n[sma$Region_n==1]="Noreste"
sma$Region_n[sma$Region_n==2]="Norte-Centro"
sma$Region_n[sma$Region_n==3]="Sur"
sma$Region_n[sma$Region_n==4]="Oeste"
sma$Region_n=factor(sma$Region_n)



set.seed(198)
samp=sample(1:141, size= 61  ,replace = F)
set.seed(254)
samp2=sample(c(1:141)[-samp], size= 61  ,replace = F)
training=sma[samp,]
validation=sma[-c(samp,samp2),]
estimation=sma[samp2,]
```


## 2. Seleccione un conjunto de modelos de regresión lineal basado en el "training set". Las variables que no se considerarán en la búsqueda serán la población total y la región. Éstas son variables protegidas por ser de interés en el análisis. Utilice distintas técnicas de selección en procesos de búsqueda automática: Si distintos criterios seleccionan distintos modelos, escoja todos los "mejores modelos".

### 2.1 Estime el modelo completo (modelo aditivo de primer orden) y calcule: Adjusted R-squared, Mallow's CP, AIC y BIC.

```{r}
modelo0=lm(TDelitos~.-PobT-Delitos-Region,data=training)
summary(modelo0)
```

El $R^{2}$ ajustado:

```{r}
summary(modelo0)$adj.r.squared
```

El criterio Mallow's Cp:

```{r}
sigma2est=sigma(modelo0)^2
ssep=sigma(modelo0)^2*summary(modelo0)$df[2]

(cp2=ssep/sigma2est-(nrow(training)-2*length(coef(modelo0))))
```
Que es el número de parametros, como era de esperarse.

El AIC:

```{r}
AIC(modelo0)
```

Y el BIC:

```{r}
BIC(modelo0)
```



### 2.2 De la librería leaps, utilice la función leaps para comparar modelos usando los métodos Adjusted R-squared y Mallow's CP. Escoja el mejor modelo en base a cada uno de los criterios.

```{r}
suppressWarnings(library(leaps))
x=model.matrix(modelo0)[,-1]
y=training[,12]

modeloleapsCP=leaps(x,y,method="Cp",nbest=1)
modeloleapsCP
bestcp=function(which,cpvec,names){
  #esta función está diseñada para cuando se usa nbest=1 en la función leaps.
  pcom=seq(1:length(cpvec))+1
  pcomrest=pcom[-length(pcom)]
  cpvecrest=cpvec[-length(cpvec)]
  if(all((cpvecrest<pcomrest)==FALSE)){
    print("Todos los Cp son mayores al número de variables. Se devuelve el modelo completo.")
    return(names[which[length(cpvec),]])
  }
  pvalid=ifelse(cpvecrest>=pcomrest,NA,cpvecrest)
  i=which.min(pcomrest-pvalid)
  if((pcomrest-pvalid)[i]>=1){
    print("Precaución: El valor del Cp está algo alejado del número de parámetros.")
  }
  cat("Índice del valor de cp en el vector de cp más cercano al número de parámetros: ",i,"\n")
  print("Variables asociadas:")
  return(names[which[i,]])
}

(vexp=bestcp(modeloleapsCP$which,modeloleapsCP$Cp,colnames(x)))

if (any(grepl("^Region_n", vexp))) {
  vexp=unique(ifelse(grepl("^Region_n", vexp), "Region_n", vexp))
}

form=as.formula(paste("TDelitos ~", paste(vexp, collapse = "+")))

xper=lm(form,data=training)

modeloleapsAR=leaps(x,y,method="adjr2",nbest=1)
(vexp2=colnames(x)[modeloleapsAR$which[which.max(modeloleapsAR$adjr2),]])

if (any(grepl("^Region_n", vexp2))) {
  vexp2=unique(ifelse(grepl("^Region_n", vexp2), "Region_n", vexp2))
}

form2=as.formula(paste("TDelitos ~", paste(vexp2, collapse = "+")))

xper2=lm(form2,data=training)
```

### 2.3 De la librería MASS utilice la función stepAIC para realizar una stepwise regression (mezclando hacia delante y hacia atrás) basado en AIC y BIC. Con el argumento scope puede indicar el modelo mínimo y el máximo. Por ejemplo para indicar que el mínimo debe incluir las variables protegidas y el máximo todas las variables puede usar scope = list(lower=~PobT+factor(Region),upper=~.). El argumento direction indica si desea realizar backward elimination, forward incluision o stepwise regression.

```{r}
library(MASS)
(xper3=stepAIC(modelo0,scope=list(lower=~Region_n,upper=~.-Region),direction="both",trace=F))

(xper4=stepAIC(modelo0,scope=list(lower=~Region_n,upper=~.-Region),direction="both",trace=F,k=log(nrow(training))))
```

### 2.4 En los modelos escogidos con los 4 métodos incluya términos de interacción y realice una nueva selección en cada modelo esta vez usando la función stepAIC con el criterio AIC.   Pista: scope = list(lower=~PobT+factor(Region),upper=~.^2).


```{r}
# newxper=update(xper, ~.^2)
# newxper2=update(xper2, ~.^2)
# newxper3=update(xper3, ~.^2)
# newxper4=update(xper4, ~.^2)

(real=stepAIC(xper,scope=list(lower=~Region_n,upper=~.^2),direction="both",trace=F))
(real2=stepAIC(xper2,scope=list(lower=~Region_n,upper=~.^2),direction="both",trace=F))
(real3=stepAIC(xper3,scope=list(lower=~Region_n,upper=~.^2),direction="both",trace=F))
(real4=stepAIC(xper4,scope=list(lower=~Region_n,upper=~.^2),direction="both",trace=F))

```


## 3. Valide todos los modelos seleccionados en la sección 2 usando el set de validación. Estime los modelos, compare los coeficientes estimados y los errores estándares con los modelos del literal 2. Calcule la media cuadrática del error de predicción y compárelo con el MSE. Escoja el modelo que tenga menor discrepancia entre MSE y MSEP.

Para el modelo basado en el criterio Mallow's Cp:

```{r}
modcp=formula(xper)
modeloCp_val=lm(modcp, data = validation)

coef(modeloCp_val)
sigma(modeloCp_val)
coef(xper)
sigma(xper)
```
Vemos que los coeficientes han cambiado, algunos bastante, como el Intercepto y PSenior, y otros menos, como Medicos y Region_nOeste.
El modelo en el set de validación tiene menor error estándar.

```{r}
H21=model.matrix(modeloCp_val)%*%solve(crossprod(model.matrix(modeloCp_val)))%*%t(model.matrix(modeloCp_val))
presspcp_val=sum((residuals(modeloCp_val)/(1-diag(H21)))^2)
msepcp_val=presspcp_val/nrow(validation)

msepcp_val-sigma(xper)^2
```

Para el modelo basado en $R^{2}$:

```{r}
modr2=formula(xper2)
modeloR2_val=lm(modr2, data = validation)
coefficients(modeloR2_val)
sigma(modeloR2_val)
coefficients(xper2)
sigma(xper2)
```
Los cambios más grandes se producen en el intercepto, y en las variables Grad y PSenior. El resto de variables mantienen bastante cercanas a las del modelo del literal 2 (aunque PCiudad ahora es negativa)

Los errores estándares son similares.

```{r}
H22=model.matrix(modeloR2_val)%*%solve(crossprod(model.matrix(modeloR2_val)))%*%t(model.matrix(modeloR2_val))
presspR2_val=sum((residuals(modeloR2_val)/(1-diag(H22)))^2)
msepR2_val=presspR2_val/nrow(validation)

msepR2_val-sigma(xper2)^2
```

Para el modelo basado en el AIC:

```{r}
modaic=formula(xper3)
modeloAIC_val=lm(modaic, data = validation)
coefficients(modeloAIC_val)
sigma(modeloAIC_val)
coefficients(xper3)
sigma(xper3)
```
Los mayores cambios están en intercepto y PSenior. Grad y PCiudad también tienen unos cuantos cambios. Los errores estándares son muy similares.

```{r}
H23=model.matrix(modeloAIC_val)%*%solve(crossprod(model.matrix(modeloAIC_val)))%*%t(model.matrix(modeloAIC_val))
presspAIC_val=sum((residuals(modeloAIC_val)/(1-diag(H23)))^2)
msepAIC_val=presspAIC_val/nrow(validation)

msepAIC_val-sigma(xper3)^2
```

Para el modelo con BIC:

```{r}
modbic=formula(xper4)
modeloBIC_val=lm(modbic, data = validation)
coefficients(modeloBIC_val)
sigma(modeloBIC_val)
coefficients(xper4)
sigma(xper4)
```

Hay ciertos cambios en el intercepto y PSenior, pero no son tan pronunciados como en modleos anteriores. El error estándar del modelo en validación es menor.

```{r}
H24=model.matrix(modeloBIC_val)%*%solve(crossprod(model.matrix(modeloBIC_val)))%*%t(model.matrix(modeloBIC_val))
presspBIC_val=sum((residuals(modeloBIC_val)/(1-diag(H24)))^2)
msepBIC_val=presspBIC_val/nrow(validation)

msepBIC_val-sigma(xper4)^2
```

La discrepancia entre la MSEP y MSE es menor en el modelo con base en el Cp.

## 4. Estime el modelo seleccionado en el literal 3 usando el set de estimación. Presente el resumen del modelo estimado.

```{r}
defi=lm(modcp,estimation)
summary(defi)
```