report2307
1.Outcome Variables
Choice of Contract & Percentage of Consumers who are still active after month T
stargazer(gp1,type='text',summary = F,title='Default Change: 09-27 ~ 10-17 vs 10-20~11-09')##
## Default Change: 09-27 ~ 10-17 vs 10-20~11-09
## ====================================================================================
## cohort num Per mon1 mon2 mon3 mon4 mon5 mon6 mon7 mon8 mon9 mon10 mon11
## ------------------------------------------------------------------------------------
## 1 before 95 0.330 0.874 0.737 0.705 0.632 0.453 0.432 0.400 0.316 0.295 0.274 0.263
## 2 before 138 0.480 1 1 0.971 0.754 0.746 0.710 0.594 0.587 0.565 0.486 0.478
## 3 before 56 0.190 1 1 1 1 1 1 1 1 1 1 0.982
## 4 after 82 0.380 0.915 0.817 0.805 0.707 0.598 0.561 0.512 0.451 0.378 0.317 0.305
## 5 after 71 0.330 1 1 1 0.789 0.789 0.732 0.676 0.648 0.648 0.521 0.465
## 6 after 63 0.290 1 1 1 1 1 1 1 0.984 0.984 0.984 0.984
## ------------------------------------------------------------------------------------stargazer(gp2,type='text',summary = F,title='Price Change: 12-20 ~ 01-09 vs 01-10~02-04')##
## Price Change: 12-20 ~ 01-09 vs 01-10~02-04
## ==================================================================
## cohort num Per mon1 mon2 mon3 mon4 mon5 mon6 mon7 mon8
## ------------------------------------------------------------------
## 1 before 82 0.440 0.829 0.744 0.634 0.573 0.512 0.463 0.439 0.354
## 2 before 40 0.220 1 1 0.875 0.825 0.800 0.775 0.700 0.700
## 3 before 63 0.340 1 0.984 0.984 0.984 0.952 0.952 0.952 0.952
## 4 after 64 0.250 0.859 0.781 0.672 0.562 0.469 0.375 0.344 0.312
## 5 after 84 0.330 1 1 0.726 0.667 0.619 0.488 0.417 0.417
## 6 after 107 0.420 1 1 1 1 0.972 0.963 0.963 0.963
## ------------------------------------------------------------------2.Engagement Score
####…specif 1
# ses.gp_trial = subset(ses.gp,trial==1)
# ses.gp_trial$med_tmr[is.na(ses.gp_trial$med_tmr)] = 0 #mean(ses.gp_trial$med_tmr,na.rm=T)
# ses.gp_trial$med_efe[is.na(ses.gp_trial$med_efe)] = 0 #mean(ses.gp_trial$med_efe,na.rm=T)
# ses.gp_trial$n_scale = scale(ses.gp_trial$n,center=F)
# ses.gp_trial$avg_problem_scale = scale(ses.gp_trial$avg_problem,center=F)
# ses.gp_trial$num_problemas_tot_scale = scale(ses.gp_trial$num_problemas_tot,center=F)
# ses.gp_trial$index = rowSums(ses.gp_trial[,c('usagerate','gamerate','med_efe','med_tmr','avg_problem_scale','n_scale','correctionrate')])
# ses.gp_trial$index = (ses.gp_trial$index - min(ses.gp_trial$index))/diff(range(ses.gp_trial$index)) *100
# ses.gp_trial$math_tarifa = factor(ses.gp_trial$math_tarifa)
# ses.gp_trial = ses.gp_trial %>% group_by(academic_year,cohort) %>% mutate(avgindex=mean(index)) %>% group_by(academic_year,cohort,math_fecha_fin) %>% mutate(avgindex_tarifa=mean(index))
#
# ses.gp_trial = ses.gp_trial %>% group_by(academic_year,cohort) %>% mutate(avgindex=mean(index)) %>% group_by(academic_year,cohort,math_tarifa) %>% mutate(avgindex_tarifa=mean(index))rcorr(as.matrix(ses.gp_trial[ses.gp_trial$academic_year=='2021-2022' & ses.gp_trial$change=='Default Change',c('usagerate','gamerate','med_efe','med_tmr','avg_problem_scale','n_scale','correctionrate')]))## usagerate gamerate med_efe med_tmr avg_problem_scale n_scale
## usagerate 1.00 0.60 0.40 -0.21 0.36 -0.21
## gamerate 0.60 1.00 0.30 -0.16 0.25 -0.09
## med_efe 0.40 0.30 1.00 0.00 0.20 0.04
## med_tmr -0.21 -0.16 0.00 1.00 -0.72 0.05
## avg_problem_scale 0.36 0.25 0.20 -0.72 1.00 0.05
## n_scale -0.21 -0.09 0.04 0.05 0.05 1.00
## correctionrate 0.78 0.52 0.30 -0.13 0.23 -0.19
## correctionrate
## usagerate 0.78
## gamerate 0.52
## med_efe 0.30
## med_tmr -0.13
## avg_problem_scale 0.23
## n_scale -0.19
## correctionrate 1.00
##
## n= 484
##
##
## P
## usagerate gamerate med_efe med_tmr avg_problem_scale n_scale
## usagerate 0.0000 0.0000 0.0000 0.0000 0.0000
## gamerate 0.0000 0.0000 0.0003 0.0000 0.0428
## med_efe 0.0000 0.0000 0.9651 0.0000 0.3360
## med_tmr 0.0000 0.0003 0.9651 0.0000 0.2586
## avg_problem_scale 0.0000 0.0000 0.0000 0.0000 0.3028
## n_scale 0.0000 0.0428 0.3360 0.2586 0.3028
## correctionrate 0.0000 0.0000 0.0000 0.0051 0.0000 0.0000
## correctionrate
## usagerate 0.0000
## gamerate 0.0000
## med_efe 0.0000
## med_tmr 0.0051
## avg_problem_scale 0.0000
## n_scale 0.0000
## correctionratechart.Correlation(ses.gp_trial[ses.gp_trial$academic_year=='2021-2022' & ses.gp_trial$change=='Default Change',c('usagerate','gamerate','med_efe','med_tmr','avg_problem_scale','n_scale','correctionrate')],histogram = T,pch=19)## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
ses.gp_trial[ses.gp_trial$academic_year=='2021-2022',] %>% ggplot(aes(x=index)) + geom_density(aes(fill=cohort),alpha=.4) + facet_wrap(~change) + geom_vline(aes(xintercept=avgindex,col=cohort)) + ggtitle('With Trial Info:Distribution of Engagement Score with Mean')+ theme_classic(base_size=10)
Kolmogorov-Smirnov Tests
a two-sample test of the null hypothesis that x and y were drawn from the same distribution.
Results:
We can not reject that the engagement/trial usage score of the ‘before’ and ‘after’ group are drawn from the same distribution for both changes.
Default Change:
we can not reject that the trial usage index for users who subscribe to 1-month/3-month contracts are drawn from the same distribution
we reject that the distribution are the same for those who subscribe to yearly contracts between before/after group
Price Change:
we can not reject that the trial usage index for users who subscribe to monthly/yearly contracts are drawn from the same distribution
we reject that the distribution are the same for those who subscribe to 3-month contracts between before/after group
ks.test(index ~ cohort, data=subset(ses.gp_trial,academic_year=='2021-2022' & change == 'Default Change'))##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: index by cohort
## D = 0.091675, p-value = 0.2711
## alternative hypothesis: two-sidedks.test(index ~ cohort, data=subset(ses.gp_trial,academic_year=='2021-2022' & change == 'Price Change'))## Warning in ks.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): p-value will be
## approximate in the presence of ties##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: index by cohort
## D = 0.10036, p-value = 0.2625
## alternative hypothesis: two-sidedses.gp_trial[ses.gp_trial$academic_year=='2021-2022' & ses.gp_trial$change=='Default Change',] %>% ggplot(aes(x=index)) + geom_vline(aes(xintercept=avgindex_tarifa,col=math_tarifa))+ geom_density(aes(fill=math_tarifa),alpha=.4)+ facet_wrap(~cohort) + ggtitle('With Trial Info:Distribution of Engagement Score')+ theme_classic(base_size=15) + coord_flip() #+ facet_wrap(~) #plot.title = element_text(size=22), 
ks.test(index ~ cohort, data=subset(ses.gp_trial,academic_year=='2021-2022' & math_tarifa==101 & change == 'Default Change'))##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: index by cohort
## D = 0.084163, p-value = 0.9016
## alternative hypothesis: two-sidedks.test(index ~ cohort, data=subset(ses.gp_trial,academic_year=='2021-2022' & math_tarifa==102 & change == 'Default Change'))##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: index by cohort
## D = 0.12238, p-value = 0.4543
## alternative hypothesis: two-sidedks.test(index ~ cohort, data=subset(ses.gp_trial,academic_year=='2021-2022' & math_tarifa==103 & change == 'Default Change'))##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: index by cohort
## D = 0.2915, p-value = 0.01072
## alternative hypothesis: two-sidedses.gp_trial[ses.gp_trial$academic_year=='2021-2022' & ses.gp_trial$change=='Price Change',] %>% ggplot(aes(x=index)) + geom_vline(aes(xintercept=avgindex_tarifa,col=math_tarifa))+ geom_density(aes(fill=math_tarifa),alpha=.4)+ facet_wrap(~cohort) + ggtitle('With Trial Info:Distribution of Engagement Score')+ theme_classic(base_size=15) + coord_flip()
ks.test(index ~ cohort, data=subset(ses.gp_trial,academic_year=='2021-2022' & math_tarifa==101 & change == 'Price Change'))##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: index by cohort
## D = 0.18401, p-value = 0.1825
## alternative hypothesis: two-sidedks.test(index ~ cohort, data=subset(ses.gp_trial,academic_year=='2021-2022' & math_tarifa==102 & change == 'Price Change'))##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: index by cohort
## D = 0.24301, p-value = 0.07691
## alternative hypothesis: two-sidedks.test(index ~ cohort, data=subset(ses.gp_trial,academic_year=='2021-2022' & math_tarifa==103 & change == 'Price Change'))##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: index by cohort
## D = 0.12516, p-value = 0.5319
## alternative hypothesis: two-sided####…specif 2
rcorr(as.matrix(ses.gp2[,c('usagerate','gamerate','med_efe','med_tmr','avg_problem_scale','n_scale','correctionrate')]))## usagerate gamerate med_efe med_tmr avg_problem_scale n_scale
## usagerate 1.00 0.64 0.86 0.64 0.77 0.07
## gamerate 0.64 1.00 0.56 0.36 0.53 0.05
## med_efe 0.86 0.56 1.00 0.72 0.67 0.04
## med_tmr 0.64 0.36 0.72 1.00 0.18 -0.01
## avg_problem_scale 0.77 0.53 0.67 0.18 1.00 0.04
## n_scale 0.07 0.05 0.04 -0.01 0.04 1.00
## correctionrate 0.81 0.57 0.69 0.52 0.60 0.01
## correctionrate
## usagerate 0.81
## gamerate 0.57
## med_efe 0.69
## med_tmr 0.52
## avg_problem_scale 0.60
## n_scale 0.01
## correctionrate 1.00
##
## n= 505
##
##
## P
## usagerate gamerate med_efe med_tmr avg_problem_scale n_scale
## usagerate 0.0000 0.0000 0.0000 0.0000 0.1437
## gamerate 0.0000 0.0000 0.0000 0.0000 0.2352
## med_efe 0.0000 0.0000 0.0000 0.0000 0.3396
## med_tmr 0.0000 0.0000 0.0000 0.0000 0.8772
## avg_problem_scale 0.0000 0.0000 0.0000 0.0000 0.3545
## n_scale 0.1437 0.2352 0.3396 0.8772 0.3545
## correctionrate 0.0000 0.0000 0.0000 0.0000 0.0000 0.8147
## correctionrate
## usagerate 0.0000
## gamerate 0.0000
## med_efe 0.0000
## med_tmr 0.0000
## avg_problem_scale 0.0000
## n_scale 0.8147
## correctionratechart.Correlation(ses.gp_trial[,c('usagerate','gamerate','med_efe','med_tmr','avg_problem_scale','n_scale','correctionrate')],histogram = T,pch=19)## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
ses.gp2 %>% ggplot(aes(x=index)) + geom_density(aes(fill=math_tarifa),alpha=.4)+ facet_wrap(~cohort) + ggtitle('With Trial Info:Distribution of Engagement Score')+ theme_classic(base_size=15) #+ coord_flip()## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
3.Choice Model
3.1.Default Change
___Multinomial Logit
Default Level: Choosing 1-month contract
Model (1) and (2) belong to the same model specification, Model (3) and (4) belong to the same model specification,…
Treatment Variable: ‘After Change’=1
week_change: Week of Subscription from the Week of Change (default level:0)
Variable ‘AvgUsage Instru’: Average Usage Rates of Individuals from the same week of subscription in the previous years to control for seasonality
stargazer(glm_0.0,glm_0.1,glm_0.2,glm_0.3,type='text',omit = c('age_cut','age_syn_ind','genero','ad'),title='Treatment: Default Change',add.lines = list(c('Control:age,gender,ad','Y','Y','Y','Y','Y','Y','Y','Y')),covariate.labels=c('After Change','week from change:-3','week from change:-2','week from change:-1','week from change:1','week from change:2','week from change:3','Trial Usage Index','AvgUsage Instru'))##
## Treatment: Default Change
## =====================================================================================================
## Dependent variable:
## -------------------------------------------------------------------------------
## 102 103 102 103 102 103 102 103
## (1) (2) (3) (4) (5) (6) (7) (8)
## -----------------------------------------------------------------------------------------------------
## After Change -0.542** 0.298 -0.539 0.605 -1.725 0.388
## (0.220) (0.252) (0.336) (0.390) (3.327) (3.508)
##
## week from change:-3 0.344 0.008 -1.531 0.452
## (0.441) (0.462) (3.763) (3.962)
##
## week from change:-2 0.520 -0.801 -1.282 -0.382
## (0.445) (0.514) (3.549) (3.752)
##
## week from change:-1 0.324 -0.975* -1.447 -0.568
## (0.452) (0.537) (3.457) (3.657)
##
## week from change:1 -0.254 -0.492 -0.227 -0.504
## (0.467) (0.496) (0.473) (0.502)
##
## week from change:2 -0.164 -0.091 -0.140 -0.107
## (0.474) (0.485) (0.480) (0.494)
##
## week from change:3 -0.024 0.215 0.084 0.171
## (0.673) (0.673) (0.743) (0.749)
##
## Trial Usage Index 3.139*** 4.717*** 3.117*** 4.525*** 3.147*** 4.704*** 3.135*** 4.515***
## (1.089) (1.373) (1.100) (1.398) (1.089) (1.371) (1.101) (1.398)
##
## AvgUsage Instru 0.286 19.709 8.154 -3.311
## (16.521) (19.075) (24.058) (25.378)
##
## Constant -1.430** -2.714*** -1.825** -2.114** -1.644 -17.136 -5.985 -0.110
## (0.700) (0.877) (0.801) (0.961) (12.114) (13.998) (14.060) (14.830)
##
## -----------------------------------------------------------------------------------------------------
## Control:age,gender,ad Y Y Y Y Y Y Y Y
## Akaike Inf. Crit. 1,029.113 1,029.113 1,037.610 1,037.610 1,031.744 1,031.744 1,041.374 1,041.374
## =====================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01___Binomial
stargazer(glm_d101.0,glm_d102.0,glm_d103.0,type='text',omit = c('age_cut','age_syn_ind','genero','ad'),dep.var.labels = c('101','102','103'))##
## ===============================================
## Dependent variable:
## -----------------------------
## 101 102 103
## (1) (2) (3)
## -----------------------------------------------
## periodafter 0.134 -0.784*** 0.905***
## (0.305) (0.295) (0.343)
##
## index -3.642*** 1.462 3.014**
## (1.008) (0.962) (1.226)
##
## avgusage_instru -6.988 -7.902 19.275
## (15.128) (14.415) (16.676)
##
## Constant 6.281 4.606 -16.763
## (11.096) (10.564) (12.235)
##
## -----------------------------------------------
## Observations 484 484 484
## Log Likelihood -297.277 -319.666 -254.042
## Akaike Inf. Crit. 614.554 659.332 528.084
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01stargazer(glm_d101.1,glm_d102.1,glm_d103.1,type='text',omit = c('age_cut','age_syn_ind','genero','ad'),dep.var.labels = c('101','102','103'),covariate.labels=c('week from change:-3','week from change:-2','week from change:-1','week from change:1','week from change:2','week from change:3','Trial Usage Index','AvgUsage Instru'))##
## =================================================
## Dependent variable:
## -----------------------------
## 101 102 103
## (1) (2) (3)
## -------------------------------------------------
## week from change:-3 -0.204 0.328 -0.192
## (0.395) (0.370) (0.387)
##
## week from change:-2 -0.094 0.850** -1.115**
## (0.405) (0.381) (0.444)
##
## week from change:-1 0.090 0.694* -1.152**
## (0.411) (0.393) (0.474)
##
## week from change:1 0.358 -0.039 -0.360
## (0.413) (0.406) (0.432)
##
## week from change:2 0.133 -0.100 -0.017
## (0.415) (0.410) (0.419)
##
## week from change:3 -0.098 -0.145 0.236
## (0.590) (0.563) (0.564)
##
## Trial Usage Index -3.553*** 1.535 2.841**
## (1.022) (0.972) (1.253)
##
## AvgUsage Instru 1.199* -1.836*** -1.823**
## (0.729) (0.711) (0.860)
##
## -------------------------------------------------
## Observations 484 484 484
## Log Likelihood -296.410 -318.017 -249.402
## Akaike Inf. Crit. 620.820 664.034 526.805
## =================================================
## Note: *p<0.1; **p<0.05; ***p<0.013.2.Price Change
___Multinomial Logit
stargazer(glm_1.0,glm_1.1,glm_1.2,glm_1.3,type='text',omit = c('age_cut','age_syn_ind','genero','ad'),title = 'Effect of Treatment: Price Change',covariate.labels=c('After Change','week from change:-3','week from change:-2','week from change:-1','week from change:1','week from change:2','week from change:3','Trial Usage Index','AvgUsage Instru'))##
## Effect of Treatment: Price Change
## ==========================================================================================
## Dependent variable:
## ----------------------------------------------------------------------
## 102 103 102 103 102 103 102 103
## (1) (2) (3) (4) (5) (6) (7) (8)
## ------------------------------------------------------------------------------------------
## After Change 0.996*** 0.701*** 0.917*** 0.629** -1.215 0.339
## (0.272) (0.243) (0.286) (0.256) (1.682) (1.478)
##
## week from change:-3 -1.105** -1.249*** -2.129 -0.903
## (0.487) (0.469) (1.495) (1.303)
##
## week from change:-2 -1.059** -0.720* -2.139 -0.374
## (0.479) (0.425) (1.558) (1.323)
##
## week from change:-1 -0.856* -0.324 -2.170 0.011
## (0.483) (0.420) (1.788) (1.567)
##
## week from change:1 0.166 0.119 0.174 0.122
## (0.436) (0.408) (0.436) (0.408)
##
## week from change:2 -0.208 -0.155 -0.324 -0.159
## (0.442) (0.410) (0.458) (0.424)
##
## week from change:3 2.262 1.746 2.244 1.684 2.258 1.732 2.227 1.682
## (1.404) (1.217) (1.414) (1.220) (1.403) (1.214) (1.416) (1.220)
##
## Trial Usage Index 7.581 6.757 10.612 0.431
## (8.884) (7.644) (11.149) (9.721)
##
## AvgUsage Instru -2.270** -1.046 -1.238 -0.278 -7.671 -5.851 -7.653 -0.926
## (0.885) (0.755) (0.918) (0.793) (6.402) (5.489) (6.408) (5.564)
##
## ------------------------------------------------------------------------------------------
## Akaike Inf. Crit. 897.761 897.761 908.550 908.550 900.698 900.698 911.473 911.473
## ==========================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01__Binomial
stargazer(glm_d101.01,glm_d102.01,glm_d103.01,type='text',omit = c('age_cut','age_syn_ind','genero','ad'),dep.var.labels = c('101','102','103'),title = 'Binomial Logit')##
## Binomial Logit
## ===============================================
## Dependent variable:
## -----------------------------
## 101 102 103
## (1) (2) (3)
## -----------------------------------------------
## periodafter -0.749*** 0.567** 0.203
## (0.232) (0.248) (0.221)
##
## index -1.912* 1.271 0.774
## (1.113) (1.226) (1.057)
##
## avgusage_instru -6.957 4.060 3.956
## (6.878) (8.043) (6.917)
##
## Constant 5.767 -5.294 -3.608
## (4.944) (5.795) (4.965)
##
## -----------------------------------------------
## Observations 417 417 417
## Log Likelihood -247.459 -234.416 -274.192
## Akaike Inf. Crit. 514.918 488.833 568.383
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01stargazer(glm_d101.11,glm_d102.11,glm_d103.11,type='text',omit = c('age_cut','age_syn_ind','genero','ad'),dep.var.labels = c('101','102','103'),covariate.labels=c('week from change:-3','week from change:-2','week from change:-1','week from change:1','week from change:2','week from change:3','Trial Usage Index','AvgUsage Instru'))##
## =================================================
## Dependent variable:
## -----------------------------
## 101 102 103
## (1) (2) (3)
## -------------------------------------------------
## week from change:-3 1.172*** -0.495 -0.756*
## (0.406) (0.429) (0.414)
##
## week from change:-2 0.851** -0.650 -0.232
## (0.388) (0.412) (0.363)
##
## week from change:-1 0.522 -0.659 0.084
## (0.389) (0.410) (0.353)
##
## week from change:1 -0.145 0.098 0.031
## (0.377) (0.347) (0.324)
##
## week from change:2 0.167 -0.111 -0.046
## (0.377) (0.360) (0.334)
##
## week from change:3 -1.888* 1.285 0.733
## (1.120) (1.235) (1.061)
##
## Trial Usage Index -0.043 -1.795** -0.502
## (0.728) (0.802) (0.690)
##
## -------------------------------------------------
## Observations 417 417 417
## Log Likelihood -246.234 -234.276 -272.272
## Akaike Inf. Crit. 518.467 494.552 570.543
## =================================================
## Note: *p<0.1; **p<0.05; ***p<0.014.Matching
library(MatchIt)
library(marginaleffects)
m.out0 = matchit(treat~ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid,data=ses.gp_con_change1,method='full',distance='glm')
summary(m.out0)##
## Call:
## matchit(formula = treat ~ ad_1 + age_syn_ind_1 + age_syn + genero_MASCULINO +
## index + num_uniq_kid + num_asso_kid, data = ses.gp_con_change1,
## method = "full", distance = "glm")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance 0.4401 0.4185 0.3015 0.9902
## ad_1 0.2222 0.1557 0.1600 .
## age_syn_ind_1 0.0324 0.0381 -0.0319 .
## age_syn 7.8009 7.7924 0.0041 0.9502
## genero_MASCULINO 0.4537 0.4602 -0.0131 .
## index 0.5389 0.5364 0.0177 0.8116
## num_uniq_kid 1.4676 1.6194 -0.2557 0.8047
## num_asso_kid 1.7500 1.9100 -0.1888 0.8782
## eCDF Mean eCDF Max
## distance 0.0846 0.1683
## ad_1 0.0665 0.0665
## age_syn_ind_1 0.0057 0.0057
## age_syn 0.0093 0.0187
## genero_MASCULINO 0.0065 0.0065
## index 0.0227 0.0819
## num_uniq_kid 0.0506 0.1024
## num_asso_kid 0.0320 0.0813
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance 0.4401 0.4401 -0.0000 1.0053
## ad_1 0.2222 0.1919 0.0728 .
## age_syn_ind_1 0.0324 0.0384 -0.0341 .
## age_syn 7.8009 7.8951 -0.0458 1.0057
## genero_MASCULINO 0.4537 0.4337 0.0401 .
## index 0.5389 0.5439 -0.0356 0.9563
## num_uniq_kid 1.4676 1.4373 0.0510 0.9791
## num_asso_kid 1.7500 1.6973 0.0622 0.9810
## eCDF Mean eCDF Max Std. Pair Dist.
## distance 0.0030 0.0231 0.0079
## ad_1 0.0303 0.0303 0.2988
## age_syn_ind_1 0.0060 0.0060 0.5305
## age_syn 0.0112 0.0369 1.1589
## genero_MASCULINO 0.0200 0.0200 0.4748
## index 0.0190 0.0611 0.4427
## num_uniq_kid 0.0133 0.0351 0.2093
## num_asso_kid 0.0169 0.0556 0.5827
##
## Sample Sizes:
## Control Treated
## All 289. 216
## Matched (ESS) 152.09 216
## Matched 289. 216
## Unmatched 0. 0
## Discarded 0. 0plot(m.out0, type = "jitter", interactive = FALSE)
plot(m.out0, type = "density", interactive = FALSE, which.xs = ~index + num_uniq_kid + age_syn)
plot(summary(m.out0))
Effect on Contract Choice
m.data0 = match.data(m.out0)
fit = glm(math_tarifa_101 ~ treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=m.data0,weights = weights, family='binomial')## Warning in eval(family$initialize): non-integer #successes in a binomial glm!avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## treat 1 - 0 0.0439 0.0451 0.974 0.33 1.6 -0.0445 0.132
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.highavg_predictions(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## treat Estimate Pr(>|z|) S 2.5 % 97.5 %
## 0 0.328 <0.001 15.1 0.259 0.406
## 1 0.372 <0.001 11.7 0.309 0.441
##
## Columns: treat, estimate, p.value, s.value, conf.low, conf.highfit = glm(math_tarifa_102 ~ treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=m.data0,weights = weights, family='binomial')## Warning in eval(family$initialize): non-integer #successes in a binomial glm!avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## treat 1 - 0 -0.153 0.0494 -3.1 0.00192 9.0 -0.25 -0.0564
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.highavg_predictions(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## treat Estimate Pr(>|z|) S 2.5 % 97.5 %
## 0 0.479 0.596 0.7 0.403 0.556
## 1 0.317 <0.001 24.8 0.261 0.379
##
## Columns: treat, estimate, p.value, s.value, conf.low, conf.highfit = glm(math_tarifa_103 ~ treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=m.data0,weights = weights, family='binomial')## Warning in eval(family$initialize): non-integer #successes in a binomial glm!avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## treat 1 - 0 0.11 0.0385 2.87 0.0041 7.9 0.035 0.186
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.highavg_predictions(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## treat Estimate Pr(>|z|) S 2.5 % 97.5 %
## 0 0.156 <0.001 49.0 0.109 0.219
## 1 0.277 <0.001 29.0 0.219 0.344
##
## Columns: treat, estimate, p.value, s.value, conf.low, conf.highfit = glm(math_tarifa_101 ~ treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=m.data0,weights = weights, family=quasibinomial())
avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights',comparison='lnratioavg',transform='exp')##
## Term Contrast Estimate Pr(>|z|) S 2.5 % 97.5 %
## treat ln(mean(1) / mean(0)) 1.13 0.335 1.6 0.881 1.45
##
## Columns: term, contrast, estimate, p.value, s.value, conf.low, conf.high, predicted, predicted_hi, predicted_lofit = glm(math_tarifa_102 ~ treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=m.data0,weights = weights, family=quasibinomial())
avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights',comparison='lnratioavg',transform='exp')##
## Term Contrast Estimate Pr(>|z|) S 2.5 % 97.5 %
## treat ln(mean(1) / mean(0)) 0.682 0.00178 9.1 0.537 0.867
##
## Columns: term, contrast, estimate, p.value, s.value, conf.low, conf.high, predicted, predicted_hi, predicted_lofit = glm(math_tarifa_103 ~ treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=m.data0,weights = weights, family=quasibinomial())
avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights',comparison='lnratioavg',transform='exp')##
## Term Contrast Estimate Pr(>|z|) S 2.5 % 97.5 %
## treat ln(mean(1) / mean(0)) 1.61 0.00716 7.1 1.14 2.28
##
## Columns: term, contrast, estimate, p.value, s.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo?.Effect on Survival Propobability
df.temp$active_dummy = ifelse(is.na(df.temp$day),0,1)
glm_act2 = glm(active_dummy~treat,data=subset(df.temp,month_contract==2),family='binomial')
glm_act3 = glm(active_dummy~treat,data=subset(df.temp,month_contract==3),family='binomial')
glm_act4 = glm(active_dummy~treat,data=subset(df.temp,month_contract==4),family='binomial')
glm_act5 = glm(active_dummy~treat,data=subset(df.temp,month_contract==5),family='binomial')
glm_act6 = glm(active_dummy~treat,data=subset(df.temp,month_contract==6),family='binomial')
glm_act7 = glm(active_dummy~treat,data=subset(df.temp,month_contract==7),family='binomial')
glm_act8 = glm(active_dummy~treat,data=subset(df.temp,month_contract==8),family='binomial')
glm_act9 = glm(active_dummy~treat,data=subset(df.temp,month_contract==9),family='binomial')
mod = list('month2'=glm_act2,'month3'=glm_act3,'month4'=glm_act4,'month5'=glm_act5,'month6' = glm_act6,'month7'= glm_act7,'month8'= glm_act7,'month9'= glm_act9)
modelsummary(mod,estimate='{estimate}{stars}',statistic = c('conf.int','p={p.value}'),gof_omit = '.*',title='Probability of Active on Month T')
| month2 | month3 | month4 | month5 | month6 | month7 | month8 | month9 |
|---|---|---|---|---|---|---|---|---|
(Intercept) | 2.905*** | 2.193*** | 1.198*** | 1.013*** | 0.793*** | 0.472*** | 0.472*** | 0.229+ |
[2.422, 3.469] | [1.828, 2.598] | [0.931, 1.478] | [0.757, 1.279] | [0.548, 1.047] | [0.237, 0.712] | [0.237, 0.712] | [-0.002, 0.463] | |
p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=0.053 | |
treat | 0.353 | 0.402 | 0.315 | 0.351 | 0.406* | 0.483* | 0.483* | 0.281 |
[-0.500, 1.279] | [-0.234, 1.076] | [-0.122, 0.762] | [-0.067, 0.777] | [0.007, 0.814] | [0.105, 0.867] | [0.105, 0.867] | [-0.078, 0.643] | |
p=0.430 | p=0.226 | p=0.162 | p=0.103 | p=0.048 | p=0.013 | p=0.013 | p=0.126 |
m.data0$active_m3 = ifelse(m.data0$day>30*3+10,1,0)
m.data0$active_m6 = ifelse(m.data0$day>30*6+10,1,0)
fit = glm(active_m3~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=m.data0,weights = weights, family='binomial')## Warning in eval(family$initialize): non-integer #successes in a binomial glm!avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## treat 1 - 0 0.0905 0.0426 2.13 0.0335 4.9 0.00707 0.174
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.highavg_predictions(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## treat Estimate Pr(>|z|) S 2.5 % 97.5 %
## 0 0.826 <0.001 48.1 0.763 0.875
## 1 0.855 <0.001 55.7 0.797 0.899
##
## Columns: treat, estimate, p.value, s.value, conf.low, conf.highfit = glm(active_m6~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=m.data0,weights = weights, family='binomial')## Warning in eval(family$initialize): non-integer #successes in a binomial glm!avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## treat 1 - 0 0.109 0.0481 2.26 0.0238 5.4 0.0145 0.203
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.highavg_predictions(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1),wts='weights')##
## treat Estimate Pr(>|z|) S 2.5 % 97.5 %
## 0 0.613 0.00649 7.3 0.532 0.688
## 1 0.727 < 0.001 31.3 0.662 0.783
##
## Columns: treat, estimate, p.value, s.value, conf.low, conf.highfit = glm(active_m6~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=subset(m.data0,math_tarifa==101),weights = weights, family='binomial')## Warning in eval(family$initialize): non-integer #successes in a binomial glm!avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1&math_tarifa==101),wts='weights')##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## treat 1 - 0 0.0365 0.083 0.439 0.66 0.6 -0.126 0.199
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.highavg_predictions(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1&math_tarifa==101),wts='weights')##
## treat Estimate Pr(>|z|) S 2.5 % 97.5 %
## 0 0.500 0.998 0.0 0.353 0.646
## 1 0.455 0.422 1.2 0.350 0.564
##
## Columns: treat, estimate, p.value, s.value, conf.low, conf.highfit = glm(active_m6~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid),data=subset(m.data0,math_tarifa==102),weights = weights, family='binomial')## Warning in eval(family$initialize): non-integer #successes in a binomial glm!avg_comparisons(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1&math_tarifa==102),wts='weights')##
## Term Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
## treat 1 - 0 0.159 0.0732 2.17 0.0302 5.0 0.0152 0.302
##
## Columns: term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.highavg_predictions(fit,variables = 'treat',vcov = ~subclass,newdata = subset(m.data0,treat==1&math_tarifa==102),wts='weights')##
## treat Estimate Pr(>|z|) S 2.5 % 97.5 %
## 0 0.611 0.0677 3.9 0.492 0.718
## 1 0.794 <0.001 14.9 0.671 0.879
##
## Columns: treat, estimate, p.value, s.value, conf.low, conf.highEffect on Monthly Usage Rates
linear regression: no covariates and no matching
lm1_1 = lm(usagerate~treat,data=subset(df.temp,month_contract==1))
lm1_2 = lm(usagerate~treat,data=subset(df.temp,month_contract==2))
lm1_3 = lm(usagerate~treat,data=subset(df.temp,month_contract==3))
lm1_4 = lm(usagerate~treat,data=subset(df.temp,month_contract==4))
lm1_5 = lm(usagerate~treat,data=subset(df.temp,month_contract==5))
lm1_6 = lm(usagerate~treat,data=subset(df.temp,month_contract==6))
lm1_7 = lm(usagerate~treat,data=subset(df.temp,month_contract==7))
lm1_8 = lm(usagerate~treat,data=subset(df.temp,month_contract==8))
lm1_9 = lm(usagerate~treat,data=subset(df.temp,month_contract==9))
modelsummary(list(lm1_1,lm1_2,lm1_3,lm1_4,lm1_5,lm1_6,lm1_7,lm1_8,lm1_9),estimate='{estimate}{stars}',stars=T,gof_omit = '.',statistic = c('conf.int','p={p.value}'),title='Effect on Monthly Usage Rates Over Time')
| (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) |
|---|---|---|---|---|---|---|---|---|---|
(Intercept) | 0.666*** | 0.554*** | 0.456*** | 0.439*** | 0.401*** | 0.376*** | 0.336*** | 0.301*** | 0.250*** |
[0.637, 0.695] | [0.517, 0.591] | [0.417, 0.494] | [0.396, 0.481] | [0.359, 0.443] | [0.334, 0.418] | [0.294, 0.377] | [0.260, 0.342] | [0.212, 0.288] | |
p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | p=<0.001 | |
treat | 0.031 | 0.011 | 0.064* | 0.015 | 0.041 | 0.027 | 0.015 | 0.007 | -0.004 |
[-0.014, 0.075] | [-0.046, 0.067] | [0.005, 0.123] | [-0.049, 0.080] | [-0.023, 0.106] | [-0.037, 0.091] | [-0.049, 0.079] | [-0.055, 0.070] | [-0.061, 0.054] | |
p=0.178 | p=0.709 | p=0.034 | p=0.639 | p=0.209 | p=0.409 | p=0.643 | p=0.821 | p=0.897 |
linear regression: with covariates and no matching
lm1_1 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==1))
lm1_2 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==2))
lm1_3 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==3))
lm1_4 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==4))
lm1_5 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==5))
lm1_6 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==6))
lm1_7 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==7))
lm1_8 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==8))
lm1_9 = lm(usagerate~treat+ad+age_syn_ind+age_syn+num_uniq_kid+num_asso_kid+genero+index+no_trialinfo,data=subset(df.temp,month_contract==9))
modelsummary(list(lm1_1,lm1_2,lm1_3,lm1_4,lm1_5,lm1_6,lm1_7,lm1_8,lm1_9),estimate='{estimate}{stars}',coef_omit = c(1,3:10),stars=T,gof_omit = '.',statistic = c('conf.int','p={p.value}'),title='Effect on Monthly Usage Rates Over Time')
| (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) |
|---|---|---|---|---|---|---|---|---|---|
treat | 0.045* | 0.023 | 0.075** | 0.027 | 0.056+ | 0.039 | 0.027 | 0.019 | 0.008 |
[0.005, 0.085] | [-0.029, 0.076] | [0.018, 0.131] | [-0.036, 0.090] | [-0.006, 0.118] | [-0.022, 0.101] | [-0.034, 0.089] | [-0.042, 0.080] | [-0.048, 0.064] | |
p=0.029 | p=0.387 | p=0.009 | p=0.396 | p=0.076 | p=0.211 | p=0.384 | p=0.542 | p=0.784 |
matching
df.temp.pvt = df.temp[,c(1:19,23,31:ncol(df.temp))] %>% pivot_wider(names_from = month_contract,values_from = usagerate,names_prefix = 'month_')
df.temp.pvt = df.temp.pvt %>% dummy_cols(select_columns = c('genero','ad','age_syn_ind'),remove_first_dummy = T)
df.temp.pvt$index[is.na(df.temp.pvt$index)] = 0
m.out1 = matchit(treat~ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo,data=df.temp.pvt,method='full',distance='glm')
summary(m.out1)##
## Call:
## matchit(formula = treat ~ ad_1 + age_syn_ind_1 + age_syn + genero_MASCULINO +
## index + num_uniq_kid + num_asso_kid + no_trialinfo, data = df.temp.pvt,
## method = "full", distance = "glm")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance 0.4292 0.4035 0.3234 1.1236
## ad_1 0.2373 0.1566 0.1898 .
## age_syn_ind_1 0.0316 0.0291 0.0146 .
## age_syn 7.9019 7.9016 0.0002 0.9672
## genero_MASCULINO 0.4684 0.4787 -0.0208 .
## index 0.5371 0.5262 0.0770 0.7604
## num_uniq_kid 1.4525 1.5839 -0.2220 0.8108
## num_asso_kid 1.7658 1.8792 -0.1319 0.9014
## no_trialinfo 0.0316 0.0582 -0.1515 .
## eCDF Mean eCDF Max
## distance 0.1087 0.2069
## ad_1 0.0807 0.0807
## age_syn_ind_1 0.0026 0.0026
## age_syn 0.0092 0.0317
## genero_MASCULINO 0.0104 0.0104
## index 0.0240 0.0612
## num_uniq_kid 0.0438 0.0880
## num_asso_kid 0.0227 0.0625
## no_trialinfo 0.0265 0.0265
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio
## distance 0.4292 0.4296 -0.0047 1.0117
## ad_1 0.2373 0.2375 -0.0004 .
## age_syn_ind_1 0.0316 0.0316 0.0003 .
## age_syn 7.9019 7.8201 0.0405 0.8775
## genero_MASCULINO 0.4684 0.4012 0.1347 .
## index 0.5371 0.5407 -0.0253 0.9621
## num_uniq_kid 1.4525 1.4556 -0.0052 0.9720
## num_asso_kid 1.7658 1.7924 -0.0309 0.9376
## no_trialinfo 0.0316 0.0393 -0.0436 .
## eCDF Mean eCDF Max Std. Pair Dist.
## distance 0.0030 0.0190 0.0126
## ad_1 0.0002 0.0002 0.1452
## age_syn_ind_1 0.0000 0.0000 0.1764
## age_syn 0.0196 0.0596 1.1369
## genero_MASCULINO 0.0672 0.0672 0.1548
## index 0.0177 0.0552 0.2856
## num_uniq_kid 0.0019 0.0044 0.2121
## num_asso_kid 0.0053 0.0127 0.1775
## no_trialinfo 0.0076 0.0076 0.1323
##
## Sample Sizes:
## Control Treated
## All 447. 316
## Matched (ESS) 175.43 316
## Matched 447. 316
## Unmatched 0. 0
## Discarded 0. 0m.data1 = match.data(m.out1)fit1 = lm(month_1~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit1_out = avg_comparisons(fit1, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit1, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
fit2 = lm(month_2~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit2_out=avg_comparisons(fit2, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit2, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
fit3 = lm(month_3~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit3_out=avg_comparisons(fit3, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit3, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
fit4 = lm(month_4~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit4_out=avg_comparisons(fit4, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit4, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
fit5 = lm(month_5~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit5_out=avg_comparisons(fit5, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit5, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
fit6 = lm(month_6~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit6_out=avg_comparisons(fit6, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit6, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
fit7 = lm(month_7~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit7_out=avg_comparisons(fit7, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit7, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
fit8 = lm(month_8~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit8_out=avg_comparisons(fit8, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit8, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
fit9 = lm(month_9~treat*(ad_1+age_syn_ind_1+age_syn+genero_MASCULINO + index + num_uniq_kid + num_asso_kid+no_trialinfo),data=m.data1,weights=weights)
fit9_out=avg_comparisons(fit9, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
#avg_predictions(fit9, variables = "treat",vcov = ~subclass,newdata = subset(m.data1, treat == 1),wts = "weights")
df = data.frame(cbind(month=1:9,rbind(fit1_out,fit2_out,fit3_out,fit4_out,fit5_out,fit6_out,fit7_out,fit8_out,fit9_out)))
df[sapply(df,is.numeric)] = lapply(df[sapply(df,is.numeric)],function(x){round(x,3)})