Corrección parcial políticas públicas
setwd("/Users/josezea/Google Drive/Laboral 2015/Docencia USTA/Políticas públicas/Parcial 2")
load("ny.Rdata")- Suponga una muestra de 10 individuos, 4 tratamientos y 6 controles, que inducen el siguiente conjunto de datos. Conformación de la base de datos:
ID <- c("A","B","C","D","E","F","G","H","I","J")
X1 <- c(41, 14, 12, 14, 13, 24, 9, 17, 5, 17)
X2 <- c(470, 139, 150, 139, 90, 90, 310, 230, 200, 230)
Z <- c(1, 0 , 0, 1, 1, 0, 0, 1, 0, 0)
Y <- c(150, 120, 180, 230, 180, 80, 75, 240, 120, 60)
punto1 <- data.frame(ID, X1, X2, Z, Y)
punto1## ID X1 X2 Z Y
## 1 A 41 470 1 150
## 2 B 14 139 0 120
## 3 C 12 150 0 180
## 4 D 14 139 1 230
## 5 E 13 90 1 180
## 6 F 24 90 0 80
## 7 G 9 310 0 75
## 8 H 17 230 1 240
## 9 I 5 200 0 120
## 10 J 17 230 0 60- Lleve a cabo una regresión logística de la variable de tratamiento Z con X1 y X2. Presente un resumen (mínimo, máximo, promedio) de las probabilidades obtenidas en el punto anterior)
m_p1 <- glm(factor(Z) ~ X1 + X2, family = "binomial", data = punto1)
summary(m_p1)##
## Call:
## glm(formula = factor(Z) ~ X1 + X2, family = "binomial", data = punto1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3685 -0.8969 -0.6507 1.1272 1.4998
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.072161 1.908372 -1.086 0.278
## X1 0.106447 0.102685 1.037 0.300
## X2 -0.000486 0.008232 -0.059 0.953
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13.460 on 9 degrees of freedom
## Residual deviance: 11.738 on 7 degrees of freedom
## AIC: 17.738
##
## Number of Fisher Scoring iterations: 4scores <- fitted.values(m_p1)
punto1$scores <- scores
summary(scores)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1629 0.3030 0.3431 0.4000 0.4075 0.8873- Construya la matriz de distancias (4 X 6) teniendo en cuenta el propensity score.
library(dplyr)##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionpunto1 <- arrange(punto1, desc(Z)) #Primero los trat y luego controles
punto1## ID X1 X2 Z Y scores
## 1 A 41 470 1 150 0.8873301
## 2 D 14 139 1 230 0.3431125
## 3 E 13 90 1 180 0.3247366
## 4 H 17 230 1 240 0.4074934
## 5 B 14 139 0 120 0.3431125
## 6 C 12 150 0 180 0.2957344
## 7 F 24 90 0 80 0.6079798
## 8 G 9 310 0 75 0.2201502
## 9 I 5 200 0 120 0.1628571
## 10 J 17 230 0 60 0.4074934D <- matrix(NA,nrow = 4, ncol = 6)
rownames(D) <- punto1$ID[1:4]
colnames(D) <- punto1$ID[5:10]
for(i in 1:4){
for(j in 1:6){
D[i,j] <- abs(punto1$scores[i] - punto1$scores[j+4] )
}
}
D## B C F G I J
## A 0.54421757 0.59159572 0.2793503 0.6671799 0.7244730 0.47983673
## D 0.00000000 0.04737815 0.2648672 0.1229623 0.1802554 0.06438084
## E 0.01837590 0.02900226 0.2832431 0.1045864 0.1618795 0.08275674
## H 0.06438084 0.11175900 0.2004864 0.1873432 0.2446362 0.00000000- Considere 0.8 desviaciones estándar como el calibrador para llevar a cabo el emparejamiento con el propensity score obtenido del punto anterior. (Considere un emparejamiento con reemplazo)
caliper <- 0.8 * sd(punto1$scores)
D<=caliper## B C F G I J
## A FALSE FALSE FALSE FALSE FALSE FALSE
## D TRUE TRUE FALSE TRUE FALSE TRUE
## E TRUE TRUE FALSE TRUE TRUE TRUE
## H TRUE TRUE FALSE FALSE FALSE TRUEEl individuo A(1) no se pudo emparejar, el individuo D(2) se emparejo con B (5), C(6) y G(8). El individuo E(3) se empreajeo con todos menos F(-7), el individuo H (4) se emparejo con B(5) y C(6)
emparejamiento <- c( c(2,5,6, 8), c(5,6,8), c(3,5,6,8,9,10), c(4,5,6))
punto1_emparejada <- punto1[emparejamiento,] - Revise el balanciamiento de X1 y X2 en la muestra emparejada, compare los resultados con la base de datos original.
Revisamos la variable X1
library(ggplot2)
ggplot(punto1, aes(factor(Z), X1)) + geom_boxplot()
ggplot(punto1_emparejada, aes(factor(Z), X1)) + geom_boxplot()
Esta variable mejora en su emparejamiento pero aún no es perfecto por lo que se sugeriría trabajar con un caliper mucho menor.
Revisamos la variable X2
ggplot(punto1_emparejada, aes(factor(Z), X2)) + geom_boxplot()
ggplot(punto1, aes(factor(Z), X2)) + geom_boxplot()
Esta variable si mejoró mucho en su emparejamiento.
- Con la base de datos completa (tabla 1) realice la regresión líneal: Y ~ Z + X1 + X2 y reporte la estimación del efecto del tratamiento, comparela este resultados con la estimación del tratamiento llevada a cabo con la base emparejada en el punto anterior.
boxplot(Y~factor(Z), data = punto1)
summary(lm(Y ~ factor(Z) + X1+ X2, data = punto1))##
## Call:
## lm(formula = Y ~ factor(Z) + X1 + X2, data = punto1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.231 -25.243 2.333 10.734 68.519
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 144.93598 30.17083 4.804 0.00299 **
## factor(Z)1 111.55058 29.46508 3.786 0.00912 **
## X1 -1.75555 1.82204 -0.964 0.37250
## X2 -0.08259 0.14634 -0.564 0.59296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 41.77 on 6 degrees of freedom
## Multiple R-squared: 0.7117, Adjusted R-squared: 0.5675
## F-statistic: 4.937 on 3 and 6 DF, p-value: 0.04638El efecto del tratamiento es 111.55058
boxplot(Y~factor(Z), data = punto1_emparejada)
summary(lm(Y ~ factor(Z) + X1+ X2, data = punto1_emparejada))##
## Call:
## lm(formula = Y ~ factor(Z) + X1 + X2, data = punto1_emparejada)
##
## Residuals:
## Min 1Q Median 3Q Max
## -61.41 -20.28 -10.05 41.17 53.95
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 207.8210 60.9492 3.410 0.00518 **
## factor(Z)1 82.2614 26.4825 3.106 0.00908 **
## X1 -1.2543 3.6471 -0.344 0.73687
## X2 -0.3596 0.1523 -2.361 0.03601 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.16 on 12 degrees of freedom
## Multiple R-squared: 0.6294, Adjusted R-squared: 0.5368
## F-statistic: 6.794 on 3 and 12 DF, p-value: 0.006269El efecto del tratamiento es 82.2614
- En los Estados Unidos existe un programa llamado colegios de excelencia que busca descubrir los colegios con metodologías de enseñana adecuadas, buenas prácticas educativas y excelentes resultados académicos para así replicarlos en los otros colegios del país. Utilize la base de datos ny.rdata en donde se presentan 584 colegios de Nueva York, las variables de la base de datos son las siguientes:
- School: nombre del colegio.
- stw: variable de agrupación; representa si el colegio ha sido asignado al programa “Colegios de excelencia”. Cada año se identifican colegios que son excelentes académicamente y en aspectos como de estructura administrativa, responsabilidad social.
- tot: tamaño del colegio
- min: porcentaje de estudiantes minoritarios en el colegio
- dis: porcentaje de estudiantes que reciben descuento en el almuerzo o lo reciben gratis.
SOLUCIÓN: puede ver este documento: http://pareonline.net/getvn.asp?v=19&n=18 Primero emparejamos con el vecino más cercano (1 control por cada tratamiento)
library(MatchIt)## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## selectmatching_p3a <- matchit(stw ~ tot + min + dis, data = ny, method = "nearest", ratio = 1)
baseemparejada_p3a <- match.data(matching_p3a) Luego emparejamos con el método de subclasificación (3 grupos) b. Realice un matching utilizando el método de subclasificación (cortes con percentiles)
m_p3<-glm(stw ~ tot + min + dis, data = ny, family="binomial")
ny$score <- m_p3$fitted.values
library(classInt)
#intervalos <-classIntervals(filter(ny,stw == 1)$score,3,style="quantile")
intervalos <- quantile(x = filter(ny,stw == 1)$score, probs = seq(0, 1, 1/3) ) #Alternativamente
names(intervalos)## [1] "0%" "33.33333%" "66.66667%" "100%"cortes <- as.numeric(intervalos)
ny$grupo <- ifelse(ny$score > cortes[1] & ny$score < cortes[2], 1,
ifelse(ny$score > cortes[2] & ny$score < cortes[3], 2,
ifelse(ny$score > cortes[3] & ny$score < cortes[4], 3, NA)))
table(ny$grupo,useNA="always")##
## 1 2 3 <NA>
## 232 54 43 255baseemparejada_p3b <- ny[!is.na(ny$grupo),]
ny$grupo <- NULL- Verifique el balanceamiento de la covariables, ¿con cuál de los dos métodos se quedaría para este problema? Revisamos los balaceamientos para las diferentes variables, empezamos con tot:
p1 <- ggplot(ny,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() +
labs(x = "Trat", y = "tot",title = "Datos sin emparejar")
p2 <- ggplot(baseemparejada_p3a,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() + labs(x = "Trat", y = "tot",title = "Vecino más cercano")
p3 <- ggplot(baseemparejada_p3b,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() + labs(x = "Trat", y = "tot",title = "Subclase (3 grupos)")
library(gridExtra)## Loading required package: gridgrid.arrange(p1, p2, p3, ncol=3)
Balanceamiento para la variable min
p4 <- ggplot(ny,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() +
labs(x = "Trat", y = "min",title = "Datos sin emparejar")
p5 <- ggplot(baseemparejada_p3a,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() + labs(x = "min", y = "tot",title = "Vecino más cercano")
p6 <- ggplot(baseemparejada_p3b,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() + labs(x = "min", y = "tot",title = "Subclase (3 grupos)")
grid.arrange(p4, p5, p6, ncol=3)
Balanceamiento para la variable min
Balanceamiento para la variable dis
p10 <- ggplot(ny,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() +
labs(x = "Trat", y = "dis",title = "Datos sin emparejar")
p11 <- ggplot(baseemparejada_p3a,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() + labs(x = "dis", y = "tot",title = "Vecino más cercano")
p12 <- ggplot(baseemparejada_p3b,aes(x =factor(stw), y = tot,fill=factor(stw))) + geom_boxplot() + geom_point() + labs(x = "dis", y = "tot",title = "Subclase (3 grupos)")
grid.arrange(p10, p11, p12, ncol=3)