taller 5 modelos

---
title: "Taller 5_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 7: Regresión Lineal Múltiple con predictores categóricos


## 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", row.names=1, quote="\"", comment.char="")
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)
```

### Codificación de las variables categóricas: Indicativas o dummy

La variable región es una variable categórica. Sin embargo, si pedimos la clase de la variable, nos saldrá como numérica:enteros.

```{r}
class(sma$Region)

```
Si se ingresa al modelo así, R va a usar a la variable región como numérica y por lo tanto va a estimar un solo coeficiente.

```{r}
summary(lm(TDelitos~Region, data=sma))

```
Para indicar que es una variable categórica vamos a utilizar la función factor().

```{r}
sma$Regionf=factor(sma$Region)
class(sma$Regionf)

#Podemos saber cúales son los niveles o clases de la variable usando levels()
levels(sma$Regionf)

```
Ahora sí, si la ingresamos en el modelo, va a esta tres coeficientes, usando una región como categoría de referencia.

```{r}
summary(lm(TDelitos~Regionf, data=sma))

```
Obtengo los mismos resultados si utilizo la función factor dentro de la fórmula del modelo.

```{r}
summary(lm(TDelitos~factor(Region), data=sma))

```

En este caso, la región 1: Noreste es la categoría de referencia.
El default de R para la categoría de referencia es usar el primer nivel (ordenado de menor a mayor, o en orden alfabético).

Por ejemplo, creemos otra variable Región 2, que en vez de usar número para las regiones usa los nombres de la región:


```{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"

class(sma$Region_n)
#La clase de la variable ahora es "character"
#Aunque se puede usar como tal, es mejor siempre utilizar factor()

sma$Region_n<-factor(sma$Region_n)

levels(sma$Region_n)

```
Las categorías están ordenados alfabéticamente, y la referencia sería Noreste

```{r}
summary(lm(TDelitos~Region_n, data=sma))

```
La interpretación de los coeficientes sería en comparación a la categoría de referencia, en este caso Noreste.

Supongamos que estamos interesados en comparar con respecto al Sur. En ese caso podemos cambiar la categoría base con la función relevel().

```{r}
relevel(x=sma$Regionf, ref = 2)

#OJO: la función relevel afecta directamente a la variable, cada vez que desee usar otra categoría de referencia debe cambiarlo con la opción ref=

```
Ahora nos sale Levels: 2 1 3 4. Lo que significa que el dos es la categoría base.


```{r}
sma$Region2=relevel(x=sma$Regionf, ref = 2)

summary(lm(TDelitos~Region2, data=sma))
```
Los coeficientes estimados cambian porque ahora son relativos a la categoría de referencia que es el Sur.

Además, podemos conocer las variables indicadoras creadas en este modelo usando contrasts()

```{r}
contrasts(sma$Region2)

```
Las columnas identifican a las tres variables indicativas. Por ejemplo, la primera variable es 1 si la región es 1 y 0 caso contrario. La primera fila que toma el valor de 0 en las tres variables es la región de referencia, en este caso la región 2.

### Pregunta 1. Estime el siguiente modelo e interprete todos los coeficientes estimados.   Interprete las pruebas t de los coeficientes de las regiones.

$$
TDelitos = \beta_0 + \beta_1 Area + \beta_2 Grad + \beta_3 Region_{Noreste} + \beta_4 Region_{Norte-Centro} + \beta_5 Region_{Sur} + \epsilon
$$
Para modelar el modelo descrito, se necesita usar como categoría de referencia a la región Oeste. Así, el modelo descrito puede ser modelado por:

```{r}
sma$Region_n=relevel(x=sma$Region_n,ref="Oeste")

modelo=lm(TDelitos~Area+Grad+Region_n,data=sma)

(res=summary(modelo))
```
$\beta_{0}=40.41$ es la tasa de crímenes promedio de una zona con 0 de área, 0 graduados, en la región Oeste.
$\beta_{1}=0.0007$ es el cambio en la tasa de crímenes promedio por cada milla cuadrada adicional, suponiendo que el porcentaje de graduados es constante, en una región cualquiera.
$\beta_{2}=0.4139$ es el cambio en la tasa de crímenes promedio por cada 1% de adicional de graduados, suponiendo que el área de la zona se mantiene igual, en una región cualquiera.
$\beta_{3}=-20.75$ es el cambio en la tasa de crímenes promedio entre la región Oeste y la región Noreste, si ambas tienen la misma área, y el mismo porcentaje de graduados.
$\beta_{4}=-12.43$ es el cambio en la tasa de crímenes promedio entre la región Oeste y la región Norte-Centro, si ambas tienen la misma área, y el mismo porcentaje de graduados.

Antes de hablar sobre la interpretación de $\beta_{5}$, discutamos la interpretación de las pruebas t. Usando un nivel de significancia del 5%, vemos que para $\beta_{i}, i=0,\cdots,4$ existe evidencia estadística para rechazar la hipótesis de que cada una de ellos es igual a 0.

Esto no ocurre para $\beta_{5}$, pues el valor de su prueba t es superior al de nuestro nivel de significancia. Por lo tanto, no existe evidencia para rechazar $H_{0}:\beta_{5}=0$. Por este motivo, no se realiza una interpretación de $\beta_{5}$.

### Pregunta 2. Ahora estime el siguiente modelo y compare la prueba F general de la regresión y el coeficiente de determinación con el modelo de la pregunta 1. ¿Existe algún cambio?

$$
TDelitos = \beta_0 + \beta_1 Area + \beta_2 Grad + \beta_3 Region_{Norte-Centro} + \beta_4 Region_{Sur} + \beta_5 Region_{Oeste} + \epsilon
$$
```{r}
sma$Region_n=relevel(x=sma$Region_n,ref="Noreste")

modelo2=lm(TDelitos~Area+Grad+Region_n,data=sma)

(res2=summary(modelo2))

res$r.squared
res2$r.squared

res$df
res$fstatistic
res2$df
res2$fstatistic
```
Los valores del $R^{2}$ y del estadístico F, así como sus grados de libertad, y por tanto, su valor p, no cambian a pesar de que usó una codifiación distinta.


### Codificación de las variables categóricas: Contraste

Para realizar la codificación por contraste, donde la referencia es el promedio de todas las categorías, podemos usar contrast() y la opción contr.sum


```{r}

contrasts(sma$Regionf)=contr.sum

contrasts(sma$Regionf)


#OJO: después de usar la opcion contr.sum, el cambio se hace directamete en la variable, por lo que desea volver a la codificación indicadora, dummy debe usar contr.treatment

```
En este caso el intercepto tendría la interpretación usando el promedio de las regiones y el coeficiente de la región 4 se obtiene como la suma de los negativos de las demás regiones.


### Pregunta 3. Modele la tasa de delitos con las variables Área, número de graduados y región. Utilice la codificación de contraste para la variable región. Interprete todos los coeficientes estimados.    Interprete las pruebas t de los coeficientes de las regiones. Compare la prueba F general de la regresión y el coeficiente de determinación con los modelos de las pregunta 1 y 2. ¿Existe algún cambio?

```{r}
contrasts(sma$Region_n)=contr.sum
contrasts(sma$Region_n)

modelo3=lm(TDelitos~Area+Grad+Region_n,data=sma)
(res3=summary(modelo3))
```

$\beta_{0}=30.74$ es la tasa de delitos promedio en una zona con 0 de área, 0 graduados, en una región promedio (esto es, sin importar si es Noereste, Norte-centro, Sur u Oeste)
$\beta_{1}=0.0007496$ es el cambio en la tasa de crímenes promedio por cada milla cuadrada adicional, suponiendo que el porcentaje de graduados es constante, en una región cualquiera.
$\beta_{2}=0.4139$ es el cambio en la tasa de crímenes promedio por cada 1% de adicional de graduados, suponiendo que el área de la zona se mantiene igual, en una región cualquiera.
$\beta_{3}=-11.08$ es el cambio en la tasa de crímenes promedio si la región es el Noreste.
$\beta_{4}=9.67$ es el cambio en la tasa de crímenes si la región es el Oeste.

Con una significancia del 5%, para los cuatro primeros $\beta$ se rechaza la hipótesis nula que indica que son iguaes a 0. Para $\beta_{5}$, dado nuestro nivel de significancia, no exeiste evidencia estadística para rechazar dicha hipótesis, por lo tanto, no se lo interpreta.

```{r}
res2$r.squared
res3$r.squared


res2$df
res2$fstatistic

res3$df
res3$fstatistic
```

De nuevo, no existe ningún cambio. Recuérdese que la prueba F y el R cuadrado del modelo de la pregunta 2 son iguales a los del de la pregunta 1.

## Más codificaciones:

Usted puede manualmente crear una codificación de interés.
Además puede crear constrastes ortogonales usando contr.poly(), contrastes de Helmert. Revise la ayuda de R ?contr.helmert. Ahí aparece una lista de todos los contrastes implementados.