DID glm

example DID glm.R

##########################################################
##########################################################
#########           under construction       #############
##########################################################
##########################################################

# see: http://econometricsense.blogspot.com/2019/01/modeling-claims-with-linear-vs-non.html?_sm_au_=iVV3RHk1r75TjsfM

# this is also the approach taken in Jonk: Yvonne Jonk, Karen Lawson, Heidi O'Connor, Kirsten S. Riise, David 
# Eisenberg, Bryan Dowd, Mary J. KreitzerMed Care. 2015 Feb; 53(2): 133–140. doi: 10.1097/MLR.0000000000000287


### OLS

summary(lm(Admit ~ treat + time + treat*time,data = did.match)) 

# Coefficients:
#  Estimate Std. Error t value            Pr(>|t|)    
# (Intercept)   1.5842     0.1760    9.00 <0.0000000000000002 ***
#  treat        -0.0693     0.2489   -0.28               0.781    
#  time         -0.5644     0.2489   -2.27               0.024 *  
#  treat:time   -0.1782     0.3520   -0.51               0.613   

### GLM
summary(m1 <- glm.nb(Admit ~ treat + time + treat*time, control = glm.control(maxit = 500),data = did.match))


# Coefficients:
#  Estimate Std. Error z value Pr(>|z|)    
#  (Intercept)   0.4601     0.1142    4.03 0.000056 ***
#  treat        -0.0447     0.1624   -0.28     0.78    
#  time         -0.4404     0.1719   -2.56     0.01 *  
#  treat:time   -0.2333     0.2500   -0.93     0.35  

# DID on model predicted values
treat <- c(1,0,1,0)
time <- c(0,0,1,1)
tmp.scoredat <- data.frame(treat,time)
tmp.scoredat$phat <- predict(m1, tmp.scoredat, type= 'response')
print(tmp.scoredat)

# treat time   phat
# 1     1    0 1.5149
# 2     0    0 1.5842
# 3     1    1 0.7723
# 4     0    1 1.0198

# implied DID based on predicted values
(.7723-1.5149)-(1.0198-1.5842) # = -.1782

# pull out predicted treatment and control pre and post values
trt_post <- tmp.scoredat[tmp.scoredat$treat == 1 & tmp.scoredat$time == 1,]
trt_pre <- tmp.scoredat[tmp.scoredat$treat == 1 & tmp.scoredat$time == 0,]
ctrl_post <- tmp.scoredat[tmp.scoredat$treat == 0 & tmp.scoredat$time == 1,]
ctrl_pre <- tmp.scoredat[tmp.scoredat$treat == 0 & tmp.scoredat$time == 0,]

# calculate DID:
(trt_post$phat - trt_pre$phat) - (ctrl_post$phat - ctrl_pre$phat)

rm(m1,treat,time,tmp.scoredat,trt_post,trt_pre,ctrl_post,ctrl_pre) # cleanup