---
title: "Cómo usar RMarkdown"
author: "Diana Sánchez"
date: "`r Sys.Date()`"
output: html_document
link-citations: TRUE
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Pregunta de investigación

¿Cuáles han sido los temas más salientes en la agenda de los movimientos
sociales?



------------------------------------------------------------------------

```{r echo=FALSE}
#install.packages("devtools")
library(devtools)
#devtools::install_github("https://github.com/ccolonescu/PoEdata")
library(PoEdata)
```

# Carga de datos

Los datos ENCUCI se pueden consultar en <htpps://inegi.org.mx>

```{r}
library(readr)
encuci <- read_csv("C:/Users/NUCSALAE/Downloads/conjunto_de_datos_ENCUCI2020_csv/conjunto_de_datos_ENCUCI_2020_VIV/conjunto_de_datos/conjunto_de_datos_ENCUCI_2020_VIV.csv")
```


```{r}
data("encuci")
names(encuci)
```

```{r}
head(cps_small)
```

```{r}
plot(cps_small$educ, cps_small$wage,
     xlab = "educación", ylab = "Salario")
```

```{r}
library(PoEdata)
data(food)
names(food)
```

```{r}
plot(food$income, food$food_exp,
     xlab = "ingreso semanal en $100", 
     ylab = "gasto semanal en alimentación en $")
```

correr un primer modelo econométrico \#### Pruebas de hipótesis

```{r}
mod1 <- lm( food$food_exp~food$income)
mod1 <- lm(food_exp~income, data=food)
attach(food)
income
mod1 <- lm(food_exp~income)
b0 <- coef(mod1)[[1]]
b1 <- coef(mod1)[[2]]
```

```{r}
plot(food$income, food$food_exp,
     xlab = "ingreso semanal en $100", 
     ylab = "gasto semanal en alimentación en $")
abline(b0,b1, col="red")
```

```{r}
summary(mod1)
```

Predicciones partimos de una muestra y buscamos hacer inferencia hacia
toda la población

```{r}
summary(income)
newx <- data.frame(income=c(10, 20, 40))
yhat <- predict(mod1, newx)
names(yhat)<- c("ingreso= $1000", "$2000", "$4000")
yhat
```

simulacion montecarlo para probar insesgadez de los coeficientes

```{r}
N <- nrow(food)
C <- 50 # numero de submuestras 
S <- 38 # el tamaño deseado de la submuestra
sumb1 <- 0
for (i in 1:C){
  set.seed(3*i)
  subsample <- food[sample(1:N, size=S, replace = TRUE),]
  mod <- lm(food_exp~income, data = subsample)
  sumb1 <- sumb1+coef(mod)[[2]]
}
print(sumb1/C, digits = 3)
```

#### Diagnóstico, ajuste y modelado

Relaciones no lineales

```{r}
data(br)
plot(br$sqft, br$price)
```

```{r}
mod <- lm(price~sqft, data=br)
b0 <- coef(mod)[[1]]
b1 <- coef(mod)[[2]]
plot(br$sqft, br$price)
abline(b0,b1, col="red")
```

```{r}
summary(mod)
```

```{r}
mod <- lm(price~I(sqft^2), data=br)
b0 <- coef(mod)[[1]]
b1 <- coef(mod)[[2]]
plot(br$sqft, br$price, col="grey")
curve(b0+b1*x^2, col="red", add = TRUE)
```

```{r}
summary(mod)
```

### Formas funcionales

#### log-lin

```{r}
library(PoEdata)
data("wa_wheat")
plot(wa_wheat$time, wa_wheat$greenough)
```

```{r}
m1 <- lm(greenough~I(time^3), data = wa_wheat)
summary(m1)
m2 <- lm(log(greenough)~time, data = wa_wheat)
summary(m2)
```

```{r}
plot(wa_wheat$time, log(wa_wheat$greenough))
```

```{r}
data("cps4_small", package = "PoEdata")
plot(cps4_small$educ, cps4_small$wage)
```

```{r}
plot(cps4_small$educ, log(cps4_small$wage))

```

#### lin-log

#### log-log

```{r}
data("newbroiler")
??newbroiler
plot(newbroiler$p, newbroiler$q)

```

```{r}
plot(newbroiler$p, log(newbroiler$q)) #log-lin
plot(log(newbroiler$p), newbroiler$q) #lin-log
plot(log(newbroiler$p), log(newbroiler$q)) #log-log Elasticidades
mod <- lm(log(q)~log(p), data = newbroiler)
summary(mod)
```

#### otro ejemplo

```{r}
data(andy)
??andy

plot(andy$price, andy$sales)
plot(andy$advert, andy$sales)
```

```{r}
#install.packages("effects")
library(effects)
mod <- lm(sales~price+advert, data = andy)
summary(mod)
```

```{r}
effprice <- effect("price", mod, confidence.level=0.9)
??effect
effprice
plot(effprice)
```

```{r}
mod <- lm(sales~price+advert+I(advert^2), data = andy)
effects <- allEffects(mod)
plot(effects)
```

#### Violaciones

##### omisión de variable relevante

```{r}
data("edu_inc", package = "PoEdata")
??edu_inc
mod1 <- lm(faminc~he+we, data=edu_inc)
summary(mod1)
```

```{r}
mod2 <- lm(faminc~he, data=edu_inc)
summary(mod2)
mod3 <- lm(faminc~he+we+kl6, data=edu_inc)
summary(mod3)
```

##### inclusión de variable irrelevante

```{r}
mod4 <- lm(faminc~he+we+kl6+xtra_x5+xtra_x6, data=edu_inc)
summary(mod4)

```

##### criterios de selección

$$AIC = ln\left(\frac{SSE}{N}\right)+\frac{2K}{N} $$

$$BIC = ln\left(\frac{SSE}{N}\right)+\frac{Kln(N)}{N} $$

```{r}
#install.packages("broom")
library(broom)
r1 <- as.numeric(glance(mod1))
r2 <- as.numeric(glance(mod2))
r3 <- as.numeric(glance(mod3))
r4 <- as.numeric(glance(mod4))
tabla <- data.frame(rbind(r1, r2, r3, r4))[,c(1,2,8,9)]
row.names(tabla)<- c("mod1", "mod2", "mod3", "mod4")
colnames(tabla) <- c("Rsq", "AdjRsq", "AIC", "BIC")
tabla
```

##### multicolinealidad

La relación intrínseca entre las variables independientes. Todo está
relacionado con todo lo demás, pero algunas cosas muestran una relación
mucho mas fuerte.

DOs tipos de multicolinealidad: 1.- Multicolinealidad perfecta (Podemos
expresar algunas variables explicativas como una combinación lineal de
las otras de manera exacta) Near Singular Matrix

2.- Multicolinealidad imperfecta (Podemos expresar algunas variables
explicativas como una combinación lineal de las otras con una variación
aleatoria)

```{r}
data("cars", package = "PoEdata")
names(cars)
mod <- lm(mpg~cyl, data = cars)
summary(mod)
mod2 <- lm(mpg~cyl+eng+wgt, data = cars)
summary(mod2)
```

```{r}
cor(cars)
```

```{r}
plot(cars$cyl, cars$eng)
```

```{r}
plot(cars$cyl, cars$wgt)
```

```{r}
plot(cars$wgt, cars$eng)
```

$$VIF_k = \frac{1}{1-R^2_k}$$

```{r}
#install.packages("car")
library(car)
vif(mod2)
```

Aquellas variables con con un VIF mayor a 10, normalmente, generan
multicolinealidad.

```{r}
mod3 <- lm(mpg~cyl+wgt, data = cars)
summary(mod3)

```

```{r}
mod4 <- lm(mpg~cyl+eng, data = cars)
summary(mod4)
```

##### heteroscedasticidad

La varianza no es la misma entre distintos niveles de las variables
explicativas

```{r}
data("food", package = "PoEdata")
mod <- lm(food_exp~income, data = food)
plot(food$income, food$food_exp, type="p",
      xlab="ingreso", ylab="gasto en alimentación")
abline(mod, col="red")
```

Método gráfico de detección

```{r}
res <- residuals(mod)
yhat <- fitted(mod)
par(mfrow = c(1, 2))
plot(food$income, res, xlab = "ingreso", ylab="residuos")
plot(yhat, res, xlab = "y estimado", ylab="residuos")
```

Prueba de heteroscedsaticidad de Breusch-Pagan

```{r}
alfa <- 0.05
ressq <- resid(mod)^2


modres <- lm(ressq~income, data = food)
N <- nobs(modres)
library(broom)
gmodres<- glance(modres)
S <- gmodres$df+1
chisqcr <- qchisq(1-alfa, S-1)
Rsqres <- gmodres$r.squared
chisq <- N*Rsqres
pval <- 1- pchisq(chisq, S-1)
```

con un valor crítico de tablas de $\chi_{cr}^2$ de 3.8 y un nivel de
confianza del 95%, el estadístico de prueba $\chi^2$ es 7.38. Con este
valor cae en la zona de rechazo de la hipótesis nula, la cual es
homoscedasticidad.

Corrección 1, por medio de los coeficientes consistentes con errores
robustos

```{r}
summary(mod)
library("car")
cov1 <- hccm(mod, type = "hc1")
#install.packages("lmtest")
library(lmtest)
food.hc1 <- coeftest(mod, vcov=cov1)
summary(mod)
```

correción 2 por medio de mínimos cuadrados ponderados

```{r}
w <- 1/food$income
food.wls <- lm(food_exp~income, weights = w, data = food)
summary(food.wls)
```

## libreria INEGI

```{r}
# Función de descarga de la ENIGH Nueva Serie (2016-2020)
enigh_nuevaserie <- function(año = NA, datos = NA){
    
    #Generales
    enig.nueva.base = "https://www.inegi.org.mx/contenidos/programas/enigh/nc/"
    nuevaserie ="NS"
    # paso1 obtener la URL ==================
    url.base = paste0(enig.nueva.base, año, "/microdatos/")
    zipdir = tempfile()
    temp.enigh = base::tempfile()
    formato_archivo = "csv"
    
    # paso 2  descargar archivos =======
    # loops para los años
    if(año==2022 |año==2020 | año == 2018 | año ==2016){nuevaserie=tolower(nuevaserie)}else {}
    ### tabla x ======
    if( datos == "viviendas" | datos == "hogares" | datos == "poblacion" | 
        datos == "gastoshogar" | datos == "erogaciones" | datos == "gastotarjetas" | 
        datos == "ingresos" | datos == "gastopersona" | datos == "trabajos"|
        datos == "agro"| datos == "noagro"| datos == "concentradohogar"){
        
        tabla = paste0("_ns_", datos,"_")
        
        url.tabla = paste0(url.base,"enigh",año,tabla,formato_archivo,".zip")
        utils::download.file(url.tabla, temp.enigh)
        # descomprima y abra
        utils::unzip(temp.enigh, exdir=zipdir)
        data.enigh.N2 = read.csv(paste0(zipdir,"/",datos,".",formato_archivo))
        names(data.enigh.N2) = tolower(names(data.enigh.N2))
        data.enigh = data.enigh.N2
        rm(data.enigh.N2)
    }
    ####
    if (exists("data.enigh")) { return(data.enigh) } else {
        stop(message("\n La opción tecleada en el parametro datos es incorrecta")) }
}
```

```{r}
concen22 <- enigh_nuevaserie(año=2022, datos = "concentradohogar")
```

```{r}
summary(concen22$ing_cor)
summary(concen22$agrope)
summary(concen22$tam_loc)
summary(concen22$deudas)
summary(concen22$medicinas)
summary(concen22$industria)
summary(concen22$factor)

```

```{r}
plot(concen22$tam_loc)
concen22$tam_loc <- as.factor(concen22$tam_loc)
mod <-lm(ing_cor~agrope+tam_loc+deudas+medicinas+industria, data=concen22 )
summary(mod)
```

## Modelo Logit

```{r}
#install.packages("aod")
library(aod)
library(ggplot2)
data <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
head(data)
```

¿cuales son los determinantes que condicionan la admisión de los
estudiantes universitarios en la UCLA?

```{r}
summary(data)
```

```{r}
sapply(data, sd)
```

Tabla de contingencia:

```{r}
xtabs(~admit+rank, data=data)
```

```{r}
data$rank <- as.factor(data$rank)
summary(data)
```

```{r}
logit <- glm(admit~gre+gpa+rank, data=data, family = "binomial")
summary(logit)
```

```{r}
confint(logit)
```

```{r}
exp(coef(logit))
```

```{r}
newdata <- with(data, data.frame(gre=mean(gre), gpa=mean(gpa), rank=factor(1:4)))
newdata
```

```{r}
newdata$rankP <- predict(logit, newdata=newdata, type="response")
newdata
```

```{r}
newdata2 <- with(data, data.frame(
  gre=rep(seq(200, 800, length.out=100),4),
  gpa=mean(gpa), rank=factor(rep(1:4, each=100))))
newdata2
```

```{r}
newdata3 <- cbind(newdata2, 
                  predict(logit, newdata=newdata2, type="link",
                          se=TRUE))
newdata3 <- within(newdata3, {
  PredicedProb <- plogis(fit)
  LL<- plogis((fit-(1.96*se.fit)))
  UL <- plogis((fit+(1.96*se.fit)))
})
head(newdata3)
```

```{r}
ggplot(newdata3,
       aes(x=gre, y=PredicedProb))+
  geom_ribbon(aes(ymin=LL, ymax=UL, fill=rank), alpha=0.2)+
  geom_line(aes(colour=rank), size=1)
```

## Modelo multinomial

```{r}
library(foreign)
library(nnet)
library(ggplot2)
#install.packages("reshape")
library(reshape)
```

El problema a evaluar es la entrada de estudiantes de preparatoria en
universidades con licenciatura con distinto enfoque (académico, general,
vocacional). Modelar la elección con base en variables de estatus
socioeconómico y calificaciones en examen de escritura

```{r}
data <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
names(data)

```

```{r}
summary(data)
```

```{r}
# tabla de contingencia
with(data, table(ses, prog))
```

```{r}
with(data, do.call(rbind, tapply(write, prog, 
                                 function(x) c(M=mean(x),
                                               SD=sd(x)))))
```

```{r}
with(data, do.call(rbind, tapply(write, ses, 
                                 function(x) c(M=mean(x),
                                               SD=sd(x)))))
```

(tarea: lo hacen para las demas materias)

modelar el logit multinomial

```{r}
data$prog2 <- relevel(data$prog, ref="academic")
test <- multinom(prog2~ses+write, data=data)
```

```{r}
summary(test)
```

```{r}
z <- summary(test)$coefficients/summary(test)$standard.errors
z
```

```{r}
p <- (1- pnorm(abs(z),0,1))*2
p
```

$$ln\left(\frac{P(prog=general)}{P(prog=academic)}\right) = b_{10}+b_{11}(ses=2)+b_{12}(ses=3)+b_{13}write $$

$$ln\left(\frac{P(prog=vocation)}{P(prog=academic)}\right) = b_{20}+b_{21}(ses=2)+b_{22}(ses=3)+b_{23}write $$

```{r}
exp(coef(test))
```

razones de riesgo relativo

Valores estimados

```{r}
head(pp <- fitted(test))
pp[1:50,]
```

```{r}
summary(data$write)
dses <- data.frame(ses=c("low", "middle", "high"),
                   write=mean(data$write))
predict(test, newdata = dses, "probs")
```

```{r}
dses <- data.frame(ses=c("low", "middle", "high"),
                   write=67)
predict(test, newdata = dses, "probs")
```

tarea: incluir y calcular probabilidades con todas las calificaciones

```{r}
dwrite <- data.frame(ses=rep(c("low", "middle", "high"), each=41),
                     write=rep(c(30:70),3))
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, 
                                  type="probs", se=TRUE))
by(pp.write[,3:5], pp.write$ses, colMeans)
```

```{r}
lpp <- melt(pp.write, id.vars = c("ses", "write"), 
            value.name="probability")
head(lpp)
```

```{r}
ggplot(lpp, aes(x=write, y=value, colour=ses))+
  geom_line()+
  facet_grid(variable~., scales="free")
```

tarea obtener gráficas para las otras calificaciones

## Modelos de Conteo

Poisson Binomial negativa poisson o binomial negativa inflada de ceros
poisson o binomial negativa truncada

### Poisson

```{r}
library(ggplot2)
#install.packages("sandwich")
library(sandwich)
#install.packages("msm")
library(msm)
```

El problea es contar el número de premios que los estudiantes obtuvieron
en la preparatoria

```{r}
data <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
names(data)
```

```{r}
summary(data)
data <- within(data, {
  prog <- factor(prog, levels=1:3, labels=c("General", 
                                            "Academic", 
                                            "Vocational"))
})
summary(data)
```

```{r}
with(data, tapply(num_awards, prog, function(x){
  sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))
}))
```

```{r}
ggplot(data, aes(num_awards, fill=prog))+
  geom_histogram(binwidth = .5, position = "dodge")
```

```{r}
m1 <- glm(num_awards~prog+math, family = "poisson", data = data)
summary(m1)
```

```{r}
cov.m1 <- vcovHC(m1, type="HC0") # calcular errores estándar robustos
std.err <- sqrt(diag(cov.m1))
r.est <- cbind(Estimate=coef(m1), 
               "Robust SE"=std.err,
               "Pr(>|z|)"= 2*pnorm(abs(coef(m1)/std.err), lower.tail = FALSE),
               LL=coef(m1)-1.96*std.err,
               UL=coef(m1)+1.96*std.err)
r.est
```

Cuenta esperada del logaritmo el coeficiente de math es 0.07

transformarlos 1) el método delta

```{r}
s  <- deltamethod(list(~exp(x1), 
                       ~exp(x2),
                       ~exp(x3),
                       ~exp(x4)),
                       coef(m1), cov.m1)
rexp.est <- exp(r.est[,-3])
rexp.est[, "Robust SE"] <- s
rexp.est
```

```{r}
s1 <- data.frame(math=mean(data$math),
                 prog=factor(1:3, 
                             levels=1:3, 
                             labels = levels(data$prog)))
s1
```

```{r}
predict(m1, s1, type = "response", se.fit = TRUE)
```

```{r}
data$phat <- predict(m1, type = "response")
data <- data[with(data, order(prog, math)),]
ggplot(data, aes(x=math, y=phat, colour = prog))+
  geom_point(aes(y=num_awards), alpha=0.5, position=position_jitter(h=.2))+
  geom_line(size=1)+
  labs(x="Calificación de matemáticas",
       y="Número esperado de premios")
```

2)  metodo 2

$$ e^\alpha$$

```{r}
ea<-exp(coef(m1)[[1]])
```

```{r}
exp(coef(m1)[[2]])*ea
```

```{r}
exp(coef(m1)[[3]])*ea
```

```{r}
exp(coef(m1)[[4]])*ea

```

```{r}
summary(data$math)
```

```{r}
exp(coef(m1)[[4]])*ea*max(data$math)
```

## Binomial negativa

Los modelos con variable de conteo binomial negativa se usan en casos de
sobredispersión

```{r}
library(foreign)
library(ggplot2)
library(MASS)
data <- read.dta("https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta")
names(data)
```

```{r}
summary(data)
data <- within(data, {
  prog <- factor(prog, levels=1:3, 
                 labels=c("General", "Academic", "Vocational"))
  id <- factor(id)
})
```

```{r}
ggplot(data, aes(daysabs, fill = prog))+
  geom_histogram(binwidth = 1)+
  facet_grid(prog~., margins=TRUE, scales="free")
```

```{r}
with(data, tapply(daysabs, prog, function(x){
  sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))
}))
```

```{r}
m1 <- glm.nb(daysabs~math+prog, data=data)
summary(m1)
```

```{r}
newdata2 <- data.frame(
  math=rep(seq(min(data$math), max(data$math), length.out=100), 3),
  prog=factor(rep(1:3, each=100), levels=1:3, labels=levels(data$prog))
)
newdata2 <- cbind(newdata2, 
                  predict(m1, newdata2, type="link", se.fit = TRUE))
newdata2 <- within(newdata2, {
  DaysAbsent<- exp(fit)
  LL <- exp(fit-1.96*se.fit)
  UL <- exp(fit+1.96*se.fit)
})
ggplot(newdata2, aes(math, DaysAbsent))+
  geom_ribbon(aes(ymin=LL, ymax=UL, fill=prog), alpha=0.25)+
  geom_line(aes(colour=prog), size=2)
```

Por ejemplo, @barro2018crecimiento dice que ....