# Matching code for India paper (Zavaleta-Cheek et al., 2022)

#libraries used: stargazer is optional to view results in htm format
library(MatchIt)
library(dplyr)
library(stargazer)

#read data
data_india<-read.csv("data_india.csv")
colnames(data_india)

#treatment: consumption of wild foods
treat<-"cons_wildfood"

#The variables of interest:
#dds: Dietary diversity score
#mds: Minimum diet score
#seven_7: Green leafy vegetables consumed last 7days
#seven_8: Vitamin A rich fruits consumed last 7days
#hours24_7: Green leafy vegetables consumed last 24 hours
#hours24_8: Vitamin A rich fruits consumed last 24 hours
interest<-c("dds","mds","seven_7","seven_8","hours24_7","hours24_8")

#continuous covariates
covs<-c("MonthlyIncome", 
        "cropdiversity_novapr", 
        "MinDistForest"
)

#covariates to perform an exact match on
exact_covs<-c("Caste","SiteName")

#months for analysis: only June and July
months_analysis <- c("June","July")

for (i in 1:length(months_analysis)){
  month<-months_analysis[i]
  print(month)
  
  #filter for each month
  data_india_month<-dplyr::filter(data_india,data_india$Month==month)
  
  #removes one row with NA for consumption of wild food
  data_india_month_complete<-data_india_month[complete.cases(data_india_month[c(covs,treat,interest)]),]
  
  #matching: notice that the algorithm of genetic matching is random and could thus lead to small differences in results
  m.out<-matchit(as.formula(paste("cons_wildfood~",paste(covs,collapse="+")))
                 ,data = data_india_month_complete
                 ,distance = "mahalanobis"
                 ,method = "genetic"
                 ,pop.size=20000
                 ,int.seed=1
                 ,unif.seed=1
                 ,ratio=1
                 ,exact = exact_covs
                 ,replace=TRUE
  )
  
  #summary of matching
  summary(m.out,standardize=TRUE)
  
  #extract matched data including matching weights
  matchit.data_india_month<-match.data(m.out)
  
  #include treatment with other covariates
  covs_forlm<-c(covs,treat)
  
  #results for dds
  dds<-glm(as.formula(paste("dds~",paste(covs_forlm,collapse="+"))),data=matchit.data_india_month,weights=matchit.data_india_month$weights,family=quasipoisson)
  
  #Food groups:
  #7: Green leafy vegetables
  #8: Vitamin A rich fruits
  
  #results for other food groups
  hrs_7<-glm(as.formula(paste("hours24_7~",paste(covs_forlm,collapse="+"))),data=matchit.data_india_month,weights=matchit.data_india_month$weights,family = quasibinomial)
  hrs_8<-glm(as.formula(paste("hours24_8~",paste(covs_forlm,collapse="+"))),data=matchit.data_india_month,weights=matchit.data_india_month$weights,family = quasibinomial)
  days_7<-glm(as.formula(paste("seven_7~",paste(covs_forlm,collapse="+"))),data=matchit.data_india_month,weights=matchit.data_india_month$weights,family = quasibinomial)
  days_8<-glm(as.formula(paste("seven_8~",paste(covs_forlm,collapse="+"))),data=matchit.data_india_month,weights=matchit.data_india_month$weights,family = quasibinomial)
  mds<-glm(as.formula(paste("mds~",paste(covs_forlm,collapse="+"))),data=matchit.data_india_month,weights=matchit.data_india_month$weights,family = quasibinomial)
  
  stargazer(dds,
            hrs_7,
            hrs_8,
            days_7,
            days_8,
            mds,
            star.char = c("+", "*", "**", "***"),
            star.cutoffs = c(.1, .05, .01, .001),
             dep.var.labels=c(
               "Dietary diversity score",
               "24 hr green leafy vegetables",
               "24 hr vit. A rich fruits",
               "Seven days green leafy vegetables",
               "Seven days vit. A rich fruits",
               "Minimum dietary diversity"),
            covariate.labels = c("Average income","Number of different crops (Nov-Apr)","Minimum distance to forest","Consumption of wild foods"),
            type="html",title=paste0("Regression outputs: ",month), digits=2,out=paste0(month,"_results after matching.htm"))
  
  
}