#install.packages("foreign")#When using a package for the first time, you have to download it first. Then you only need to library it every time when you use it.

library(foreign)
data=read.dta("C:\\NCME 2016\\RMPW\\Riverside_first_imputed_Sept2011.dta")
attach(data)
l0 = glm(emp ~ emp_prior + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[treat==0,],family=binomial)
p0 = predict(l0,data,type="response")
xb0 = predict(l0,data)
l1 = glm(emp ~ emp_prior + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[treat==1,],family=binomial)
p1 = predict(l1,data,type="response")
xb1 = predict(l1,data)

#install.packages("MASS")
library(MASS)
data$empcat = as.factor(data$empcat)
lo0 = polr(empcat ~ emp_prior + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[treat==0,])
po0 = predict(lo0,data,type="probs")
p00 = po0[,1]
p01 = po0[,2]
p02 = po0[,3]
xbo0 = cbind(emp_prior, pqtrunc50, pqtrunc51, pqtrunc52, pqtrunc53, pqtrunc30, hispanic, pqtrunc49, nevmar)%*%c(coef(lo0))

lo1 = polr(empcat ~ emp_prior + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[treat==1,])
po1 = predict(lo1,data,type="probs")
p10 = po1[,1]
p11 = po1[,2]
p12 = po1[,3]
xbo1 = cbind(emp_prior, pqtrunc50, pqtrunc51, pqtrunc52, pqtrunc53, pqtrunc30, hispanic, pqtrunc49, nevmar)%*%c(coef(lo1))

#### Run binary logit for teen-parents as moderator
for(j in 0:1){
  for(k in 0:1){
    l = glm(emp ~ emp_prior + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[treat==j & teen_parent==k,],family=binomial)
    assign(paste0("pt",j,k),predict(l,data,type="response"))
    assign(paste0("xbt",j,k),predict(l,data))
  }
}

####### Exclusion for binary case #######
for(a in 0:1){
  assign(paste0("sd",a), sd(get(paste0("xb",a))))
  for(i in 0:1){
    for(j in 0:1){
      assign(paste0("max",a,i,j), max(get(paste0("xb",a))[treat==i & emp==j])+0.2*get(paste0("sd",a)))
      assign(paste0("min",a,i,j), min(get(paste0("xb",a))[treat==i & emp==j])-0.2*get(paste0("sd",a)))
    }
  }
}

exclude = rep(0, length(data[,1]))
for(a in 0:1){
  exclude[get(paste0("xb",a)) < max(get(paste0("min",a,0,0)), get(paste0("min",a,0,1)), get(paste0("min",a,1,0)), get(paste0("min",a,1,1)))] = 1
  exclude[get(paste0("xb",a)) > min(get(paste0("max",a,0,0)), get(paste0("max",a,0,1)), get(paste0("max",a,1,0)), get(paste0("max",a,1,1)))] = 1
}
data=cbind(data,exclude)
###############################################


####### Exclusion for 3-category case #######
for(a in 0:1){
  assign(paste0("sdo",a), sd(get(paste0("xbo",a))))
  for(i in 0:1){
    for(j in 0:2){
      assign(paste0("max",a,i,j), max(get(paste0("xbo",a))[treat==i & empcat==j])+0.2*get(paste0("sdo",a)))
      assign(paste0("min",a,i,j), min(get(paste0("xbo",a))[treat==i & empcat==j])-0.2*get(paste0("sdo",a)))
    }
  }
}

exclude3 = rep(0, length(data[,1]))
for(a in 0:1){
  exclude3[get(paste0("xbo",a)) < max(get(paste0("min",a,0,0)), get(paste0("min",a,0,1)), get(paste0("min",a,0,2)), get(paste0("min",a,1,0)), get(paste0("min",a,1,1)), get(paste0("min",a,1,2)))] = 1
  exclude3[get(paste0("xbo",a)) > min(get(paste0("max",a,0,0)), get(paste0("max",a,0,1)), get(paste0("max",a,0,2)), get(paste0("max",a,1,0)), get(paste0("max",a,1,1)), get(paste0("max",a,1,2)))] = 1
}
data=cbind(data,exclude3)
######################################################


####### Exclusion for binary case with teen-parent as moderator #########
for(c in 0:1){
  for(a in 0:1){
    assign(paste0("sdt",a,c), sd(get(paste0("xbt",a,c))[teen_parent==c]))
    for(i in 0:1){
      for(j in 0:1){
        assign(paste0("xbt",a,i,j,c,"max"), max(get(paste0("xbt",a,c))[treat==i & emp==j & teen_parent==c])+0.2*get(paste0("sdt",a,c)))
        assign(paste0("xbt",a,i,j,c,"min"), min(get(paste0("xbt",a,c))[treat==i & emp==j & teen_parent==c])-0.2*get(paste0("sdt",a,c)))
      }
    }
  }
}

excludet = rep(0, length(data[,1]))
for(c in 0:1){
  for(a in 0:1){
    excludet[teen_parent==c & get(paste0("xbt",a,c)) < max(get(paste0("xbt",a,0,0,c,"min")), get(paste0("xbt",a,0,1,c,"min")), get(paste0("xbt",a,1,0,c,"min")), get(paste0("xbt",a,1,1,c,"min")))] = 1
    excludet[teen_parent==c & get(paste0("xbt",a,c)) > min(get(paste0("xbt",a,0,0,c,"max")), get(paste0("xbt",a,0,1,c,"max")), get(paste0("xbt",a,1,0,c,"max")), get(paste0("xbt",a,1,1,c,"max")))] = 1
  }
}
data=cbind(data,excludet)
#######################################################

table(exclude,treat)
table(exclude3, treat)
table(excludet, treat)

####### Generate weights ########

###### Binary #######
rmpw=rep(1, length(data[,1]))
rmpw[treat==1 & emp==1] = as.numeric(p0/p1)[treat==1 & emp==1]
rmpw[treat==1 & emp==0] = as.numeric((1-p0)/(1-p1))[treat==1 & emp==0]

########## Binary non-parametric 3x3 ##########
c=quantile(xb1, probs = c(0.33,0.67))
h1=rep(3, length(data[,1]))
h1[xb1<=c[1]]=1
h1[xb1>c[1] & xb1<=c[2]]=2

for(j in 1:3){
  c=quantile(xb0[h1==j], probs = c(0.33,0.67))
  temp=rep(3, length(data[,1]))
  temp[xb0<=c[1]]=1
  temp[xb0>c[1] & xb0<=c[2]]=2
  assign(paste0("h0",j),temp)
}

strata=rep(NA, length(data[,1]))
strata[h1==1 & h01==1]=0
strata[h1==1 & h01==2]=1
strata[h1==1 & h01==3]=2
strata[h1==2 & h02==1]=3
strata[h1==2 & h02==2]=4
strata[h1==2 & h02==3]=5
strata[h1==3 & h03==1]=6
strata[h1==3 & h03==2]=7
strata[h1==3 & h03==3]=8

for(i in 0:1){
  for(j in 0:8){
    assign(paste0("pr",i,j),mean(emp[treat==i & strata==j & exclude==0]))
  }
}

nrmpw33 = rep(1, length(data[,1]))
for(j in 0:8){
  nrmpw33[treat==1 & strata==j & emp==0] = (1 - get(paste0("pr",0,j)))/(1 - get(paste0("pr",1,j)))
  nrmpw33[treat==1 & strata==j & emp==1] = get(paste0("pr",0,j))/get(paste0("pr",1,j))
}

########## Binary non-parametric  4x4 ##########

c=quantile(xb1, probs = c(0.25,0.5,0.75))
g1=rep(4, length(data[,1]))
g1[xb1<=c[1]]=1
g1[xb1<=c[2] & xb1>c[1]]=2
g1[xb1<=c[3] & xb1>c[2]]=3

for(j in 1:4){
  c=quantile(xb0[g1==j], probs = c(0.25,0.5,0.75))
  temp=rep(4, length(data[,1]))
  temp[xb0<=c[1]]=1
  temp[xb0>c[1] & xb0<=c[2]]=2
  temp[xb0>c[2] & xb0<=c[3]]=3
  assign(paste0("g0",j),temp)
}

strata44=rep(NA, length(data[,1]))
strata44[g1==1 & g01==1]=0
strata44[g1==1 & g01==2]=1
strata44[g1==1 & g01==3]=2
strata44[g1==1 & g01==4]=3
strata44[g1==2 & g02==1]=4
strata44[g1==2 & g02==2]=5
strata44[g1==2 & g02==3]=6
strata44[g1==2 & g02==4]=7
strata44[g1==3 & g03==1]=8
strata44[g1==3 & g03==2]=9
strata44[g1==3 & g03==3]=10
strata44[g1==3 & g03==4]=11
strata44[g1==4 & g04==1]=12
strata44[g1==4 & g04==2]=13
strata44[g1==4 & g04==3]=14
strata44[g1==4 & g04==4]=15

for(i in 0:1){
  for(j in 0:15){
    assign(paste0("pn",i,j),mean(emp[treat==i & strata44==j]))
  }
}

nrmpw44 = rep(1, length(data[,1]))
for(j in 0:15){
  nrmpw44[treat==1 & strata44==j & emp==0] = (1 - get(paste0("pn",0,j)))/(1 - get(paste0("pn",1,j)))
  nrmpw44[treat==1 & strata44==j & emp==1] = get(paste0("pn",0,j))/get(paste0("pn",1,j))
}

##### 3-category ######
rmpw3 = rep(1, length(data[,1]))
for(i in 0:2){
  rmpw3[treat==1 & get(paste0("empcat",i))==1] = (get(paste0("p",0,i))/get(paste0("p",1,i)))[treat==1 & get(paste0("empcat",i))==1] 
}

###### Binary with teen-parent moderator  ########
rmpwt = rep(1, length(data[,1]))
for(k in 0:1){
  rmpwt[treat==1 & emp==1 & teen_parent==k] = (get(paste0("pt",0,k))/get(paste0("pt",1,k)))[treat==1 & emp==1 & teen_parent==k] 
  rmpwt[treat==1 & emp==0 & teen_parent==k] = ((1-get(paste0("pt",0,k)))/(1-get(paste0("pt",1,k))))[treat==1 & emp==0 & teen_parent==k]
}
if(length(which(is.na(trunc_dep12sm2)|is.na(treat)|is.na(rmpw)))>0)
data=data[-which(is.na(trunc_dep12sm2)|is.na(treat)|is.na(rmpw)),]

#### This serves to generate a unique identifier, called "obs", for each observation. Such that duplicates will have the same identifier
obs=1:length(data[,1])
data=cbind(obs,data)

######## Generate duplicate observations, where D is indicator for duplicate   #######
D=c(rep(0, length(data[,1])),rep(1,sum(treat==1)))
data=rbind(data,data[treat==1,])
data=cbind(data, D)

####### Make sure duplicates get a weight=1  #########
rmpw[D==1]=1
nrmpw33[D==1]=1
nrmpw44[D==1]=1
rmpw3[D==1]=1
rmpwt[D==1]=1

### Generate interactions  #######
treatD=data$treat*D
teen_treatD=data$teen_treat*D
not_teen= (-1*data$teen_parent) + 1
not_teen_treat=not_teen*data$treat
not_teen_treatD=not_teen_treat*D

data=cbind(data,rmpw,nrmpw33,nrmpw44,rmpw3,rmpwt,treatD,teen_treatD,not_teen,not_teen_treat,not_teen_treatD)

####### Outcome models  #########

#Binary RMPW
l = lm(trunc_dep12sm2 ~ treat + treatD + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[data$exclude==0,], weights=rmpw)
summary(l)

#Binary Non-parametric RMPW (NRMPW)
l = lm(trunc_dep12sm2 ~ treat + treatD + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[data$exclude==0,], weights=nrmpw33)
summary(l)
l = lm(trunc_dep12sm2 ~ treat + treatD + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[data$exclude==0,], weights=nrmpw44)
summary(l)

# 3-category employment
l = lm(trunc_dep12sm2 ~ treat + treatD + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[data$exclude3==0,], weights=rmpw3)
summary(l)

# Binary RMPW with teen-parent as moderator
l = lm(trunc_dep12sm2 ~ treat + treatD + teen_parent + teen_treat + teen_treatD + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + teen_pqtrunc50 + teen_pqtrunc51 + teen_pqtrunc52 + teen_pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[data$excludet==0,], weights=rmpwt)
summary(l)
l = lm(trunc_dep12sm2 ~ treat + treatD + not_teen + not_teen_treat + not_teen_treatD + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + teen_pqtrunc50 + teen_pqtrunc51 + teen_pqtrunc52 + teen_pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[data$excludet==0,], weights=rmpwt)
summary(l)

#Binary RMPW and NRMPW with a binary depresssion measure as the outcome
l = lm(depsmrk2 ~ treat + treatD + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[data$exclude==0&!is.na(data$depsmrk2),], weights=rmpw)
summary(l)
l = lm(depsmrk2 ~ treat + treatD + pqtrunc50 + pqtrunc51 + pqtrunc52 + pqtrunc53 + pqtrunc30 + hispanic + pqtrunc49 + nevmar, data = data[data$exclude==0&!is.na(data$depsmrk2),], weights=nrmpw44)
summary(l)

# SEM/Path Analysis approach
l = lm(emp ~ treat, data = data[D==0,])
summary(l)
mediator = coef(l)["treat"]
mediator_se = coef(summary(l))[2, "Std. Error"]

l = lm(trunc_dep12sm2 ~ treat + emp, data = data[D==0,])
summary(l)
direct = coef(l)["treat"]
direct_se = coef(summary(l))[2, "Std. Error"]
indirect = coef(l)["emp"]*mediator
indirect_se = sqrt((mediator^2)*coef(summary(l))[3, "Std. Error"]^2 + (mediator_se^2)*coef(l)["emp"]^2)

direct
direct_se
indirect
indirect_se