Code: Bulk regressions

See output here.

Download R script here.

library(dplyr)
library(texreg)

#download data
t<-tempfile()
tt<-tempdir()
download.file("https://doi.org/10.1371/journal.pone.0136402.s001",destfile=t)
unzip(zipfile=t,files="blacklist_data_20150623.csv",exdir=tt)
#read
d<-read.csv(file.path(tt,"blacklist_data_20150623.csv"))
#new variables
d$id<-d$municip_geocode_id
d$treat<-I(d$emb_cont_share)*1
#minimal data set
vlist<-c(
  "id","year","def_area","treat","gdppc_l1","nibamafines_no_l1"
)
d<-d[,vlist]
#logarithmic transformations (inv. hyp. sine)
for (v in c("def_area","treat","gdppc_l1","nibamafines_no_l1")){
  d[,paste0("asinh.",v)]<-asinh(d[,v])
}


########################
### sample selection ###
########################
##all
d.all<-d
##subset
d<-d.all
i<-which(d$year>=2004)
d<-d[i,]
d.base<-d

##########################
### estimation formula ###
##########################
f.estimation<-function(m){
    # m<-modela
    ##FE normal
    if (m$model=="fe"){
      f<-as.formula(paste(m$dep,"~",paste(get(m$cov),collapse="+"),"| ",m$fe," | ",m$ivs," | ",m$cluster))
      dd<-get(m$sample)
      e<-felm(f,data=dd)
      # s<-summary(e)
    }
    ##FE weighted by deforestation probablity
    if (m$model=="fe.wp"){
      # print("fe weighted by pre-def-prob")
      # wtsin<-x$data$weight * x$data$predefprob
      # e<-felm(x$formula,data=x$data,weights=wtsin)
      # s<-summary(e)
    }
    ##RE normal
    if (m$model=="re"){
      # print("re")
      # e<-plm::plm(x$formula,data=x$data,model="random")
      # e<-plm::plm(f,data=x$data,model="random")
      # s<-summary(e)
    }

    return(e)
}


######################################
##### first step single regression ###
######################################
{
#covariate list
cov.01<-c("treat","asinh.gdppc_l1","asinh.nibamafines_no_l1")
##set parameters
m<-list(
  #dependent
  "asinh.def_area",
  #covariates
  "cov.01",
  #sample
  "d.all",
  #model
  "fe",
  #fixed effects
  "id + year",
  #ivs
  "0",
  #cluster variables
  "id"
)
names(m)<-c("dep","cov","sample","model","fe","ivs","cluster")
m
#model specification
model.first<-m

#estimate
e<-f.estimation(model.first)
summary(e)
}

######################################################
#### looping over covariate and dependent list #######
######################################################
{
##set of covariates
#delete previous covariates
ls()[grep("cov.[0-9][0-9]",ls())]
rm(list=ls()[grep("cov.[0-9][0-9]",ls())])
#set covariates
cov.01<-c("treat","asinh.gdppc_l1")
cov.02<-c("treat","asinh.gdppc_l1","asinh.nibamafines_no_l1")

##models
g<-list(
  list("asinh.def_area","cov.01","d.all","fe","id + year","0","id"),
  list("asinh.def_area","cov.02","d.all","fe","id + year","0","id"),
  list("log(def_area+1)","cov.01","d.all","fe","id + year","0","id"),
  list("log(def_area+1)","cov.02","d.all","fe","id + year","0","id")
)
#naming parameters
g<-lapply(g,function(m){names(m)<-c("dep","cov","sample","model","fe","ivs","cluster");return(m)})
#naming models (could be put in a function)
n<-unlist(lapply(g,function(x){x[["dep"]]}))
n<-gsub("\\(",".",n)
n<-gsub("\\)",".",n)
n<-gsub("\\.$","",n)
n<-gsub("\\+","",n)
n1<-n
n<-unlist(lapply(g,function(x){x[["sample"]]}))
n2<-n
n<-unlist(lapply(g,function(x){x[["cov"]]}))
n3<-n
n<-cbind(n1,n2,n3)
n<-paste(n[,1],n[,2],n[,3],sep="-")
names(g)<-n
modelgroup<-g

##estimating
elist<-lapply(
  modelgroup,
  f.estimation
)
lapply(elist,summary)

}



################################
### html output in one table ###
################################
#create table
a<-htmlreg(
  elist,
  stars=c(0.01,0.05,0.1),
  digits=3
)
html.link<-tempfile(fileext=".html")
cat(
  file=html.link,
  a
)
clipr::write_clip(html.link)
file.copy(html.link,"bulk.est.newest.html",overwrite=T)
file.copy(html.link,"bulk.est.onetable.html",overwrite=T)


##################################################
### html output in stacked tables by dependent ###
##################################################
##group models by dependent
g<-names(elist)
g<-stringr::str_split(g,"-")
g<-unlist(lapply(g,function(s){s[[1]]}))
gi<-as.integer(factor(g))
gelist<-split(elist,gi)

##rename models
names(gelist)<-unique(g)
gelist<-lapply(
  names(gelist),
  function(x){
    # x<-names(gelist)[1]
    mg<-gelist[[x]]
    n<-gsub(x,"",names(mg))
    n<-gsub("^-","",n)
    names(mg)<-n
    return(mg)
  }
)
names(gelist)<-unique(g)

#create stacked tables
a<-lapply(
  names(gelist),
  function(x){
    # x<-names(gelist)[1]
    mg<-gelist[[x]]
    html.group<-c(
      x,
      htmlreg(
        mg,
        stars=c(0.01,0.05,0.1),
        digits=3
      )
    )
    return(html.group)
  }
)
a<-unlist(a)
html.link<-tempfile(fileext=".html")
cat(
  file=html.link,
  a
)
clipr::write_clip(html.link)
file.copy(html.link,"bulk.est.newest.html",overwrite=T)
file.copy(html.link,"bulk.est.onetable.html",overwrite=T)