Northwestern GRPL Assessment
Tanvi Tilloo
2025-02-25
View the datasets and look at the column names
view(assets)
colnames(assets)## [1] "hhid" "wave" "Asset_Type" "InstanceNumber"
## [5] "animaltype" "toolcode" "durablegood_code" "quantity"
## [9] "currentvalue"view(demographics)
colnames(demographics)## [1] "wave" "hhid" "hhmid" "villageid"
## [5] "treat_hh" "gender" "age" "relationship"
## [9] "maritalstatus" "spouseinhouse" "agemarried" "religion"
## [13] "religionother" "fatherinhouse" "fathereduc" "fathereducother"
## [17] "motherinhouse" "mothereduc" "mothereducother"view(depression)
colnames(depression)## [1] "wave" "hhid" "hhmid" "tired"
## [5] "nervous" "sonervous" "hopeless" "restless"
## [9] "sorestless" "depressed" "everythingeffort" "nothingcheerup"
## [13] "worthless"Part I
Q1:
At what level is each dataset uniquely identified (i.e., what does each row represent, and which variables are unique identifiers)?
\(\color{red}{\text{SOLUTION:}}\)
The assets data is uniquely identified at the asset level and is represented by the instance number variable. Each row contains information about unique assets owned by the household in the survey. The demographics data is unique at the household member level and is given by the hhmid variable. Each row in this dataset contains information about the demographics of unique household members alongwith household id and surveyor information that ties the household member to the sampled household. Lastly, the depression data is also unique at the household member level and is represented by hhmid variable. Each row in this dataset contains a unique household member’s responses to the Kessler Psychological Distress Scale questions.
Demographics
Q2:
Import the demographics dataset and calculate a proxy variable for household size based on the number of members surveyed in each household in Wave 1. Assume the household size for Wave 2 remains the same.
#check for na values in wave
any(is.na(demographics$wave))## [1] FALSE#filter data by wave 1 and then group by household id to proxy household size
demographics_hhsize <- demographics %>% filter(wave == 1)%>% group_by(hhid) %>% summarize(hhsize = n())
#consider joining demographics_hhszie to demographics on hhid
demographics_proxy <- demographics %>% left_join(demographics_hhsize, by = "hhid")
#view(demographics_hhsize)\(\color{blue}{\text{We proxied household size by summing the number of times members from the same household id were surveyed in Wave 1 of the survey. After filtering the data to only keep observations from Wave 1, we grouped the data by household id and summed the total number of times the household id appeared in the dataset.}}\)
Assets
Q3
To calculate the monetary value of all assets, you should use the ‘currentvalue’ variable, which reports the monetary value of a single unit of the asset. However, you will notice that this variable is often missing. Please use the median of “currentvalue” for each type of asset (by type we mean, for example, “chickens”, “Cutlass”, “Room Furniture”, “Radio”, “Cell (mobile) Phone handset”, etc.) to impute the missing values.
#check for missing values in the hhid column to make sure we
#have household ids for each asset
any(is.na(assets$hhid))## [1] FALSE#check for missing values in the main asset type column to make sure we
#only work on available data
any(is.na(assets$Asset_Type))## [1] FALSE#first convert the asset type columns to categorical variables
assets <- assets %>% mutate(across(c(animaltype, toolcode, durablegood_code),
as.factor)) %>% mutate(across(c(hhid, wave), as.double))
#group by asset category and calculate the median current value for each asset type within the category
median_values <- assets %>% group_by(animaltype, toolcode, durablegood_code) %>%
summarize(median_currentvalue = median(currentvalue, na.rm = TRUE, .groups = "drop"))
#fill in missing values in current value with the corresponding median current
#value of each asset type
assets <- assets %>%
left_join(median_values, by = c("animaltype", "toolcode", "durablegood_code")) %>%
mutate(currentvalue = ifelse(is.na(currentvalue), median_currentvalue, currentvalue))Q4
Create a variable that contains the total monetary value for each observation, by multiplying quantity and the imputed current value.
assets <- assets %>% mutate(total_value = quantity * currentvalue)Q5
Produce a dataset at the household-wave level (for each household, there should be at most two observations, one for each wave) which contains the following variables:
Household ID
Wave ID
Total value of animals
Total value of tools
Total value of durable goods
Total asset value
#group by household id and wave
#then sum current value over non-na entries of asset type
assets_hhwave <- assets %>%
group_by(hhid, wave) %>%
summarize(
total_value_animals = sum(currentvalue[!is.na(animaltype)], na.rm = TRUE),
total_value_tools = sum(currentvalue[!is.na(toolcode)], na.rm = TRUE),
total_value_durable_goods = sum(currentvalue[!is.na(durablegood_code)], na.rm = TRUE),
total_asset_value = sum(currentvalue, na.rm = TRUE),
.groups = "drop"
)
#rename the columns
names(assets_hhwave) <- c("Household ID", "Wave ID", "Total value of
animals", "Total value of tools", "Total value
of durable goods", "Total asset value")Mental Health/Depression Data
Q6:
A Kessler-10 scale is a measure of mental health that uses 10 questions that identify how often people experience symptoms associated with depression. Using this reference, construct the kessler score (name it kessler score) and a categorical variable named kessler categories with 4 categories: no significant depression, mild depression, moderate depression, and severe depression. Be careful: sometimes variables in the Kessler section are missing, which can complicate construction of a score. Please justify, using comments in your code, how you handle missing values.
\(\color{red}{\text{SOLUTION:}}\) To handle missing values in the Kesler scale questions, I would examine the correlation between questions to see which responses are highly correlated. If, say, higher values on questions related to tiredness are correlated with higher values on questions related to worthlessness then we can use the available values in either of these questions to fill in missing values. Similarly, questions that follow each other in measuring the degree of a certain feeling (e.g. columns like “nervous” and “sonervous”) are also likely to be highly correlated and thus, worth exploring as an imputation strategy}}$
#first find which columns have missing values
colSums(is.na(depression))## wave hhid hhmid tired
## 0 0 0 31
## nervous sonervous hopeless restless
## 34 40 34 88
## sorestless depressed everythingeffort nothingcheerup
## 48 45 60 94
## worthless
## 88As we can see, several columns have missing values with “restless”, “nothingcheerup”, and “worthless” columns having the highest amount of missing data. These numbers are high enough that if we dropped all observations with missing values we would lose a lot of data. Thus, imputing might be necessary. At first pass, if we believe that there are correlations between a person’s responses across the survey then we can use the most correlated variables (values > 0.7) to impute missing values. Thus, we can create a correlation matrix to look for highly correlated columns.
colnames(depression)## [1] "wave" "hhid" "hhmid" "tired"
## [5] "nervous" "sonervous" "hopeless" "restless"
## [9] "sorestless" "depressed" "everythingeffort" "nothingcheerup"
## [13] "worthless"cor_matrix <- cor(depression[, c("tired", "nervous", "sonervous", "hopeless",
"restless", "sorestless", "depressed",
"everythingeffort", "nothingcheerup", "worthless")],
use = "pairwise.complete.obs")
print(cor_matrix)## tired nervous sonervous hopeless restless
## tired 1.0000000 0.44240321 0.30501116 0.29242706 0.3371423
## nervous 0.4424032 1.00000000 0.52866983 0.33226797 0.3693095
## sonervous 0.3050112 0.52866983 1.00000000 0.31173037 0.3744391
## hopeless 0.2924271 0.33226797 0.31173037 1.00000000 0.3092811
## restless 0.3371423 0.36930953 0.37443908 0.30928110 1.0000000
## sorestless 0.2859619 0.36763325 0.51182602 0.31084450 0.5589534
## depressed 0.3351592 0.44043525 0.35689863 0.41870957 0.3479190
## everythingeffort 0.1351178 0.09870054 0.04675862 0.09886671 0.1738176
## nothingcheerup 0.2981887 0.34198638 0.40501239 0.39277064 0.3391959
## worthless 0.2294393 0.30958094 0.32127135 0.50510592 0.2932832
## sorestless depressed everythingeffort nothingcheerup worthless
## tired 0.2859619 0.3351592 0.13511783 0.2981887 0.2294393
## nervous 0.3676333 0.4404353 0.09870054 0.3419864 0.3095809
## sonervous 0.5118260 0.3568986 0.04675862 0.4050124 0.3212714
## hopeless 0.3108445 0.4187096 0.09886671 0.3927706 0.5051059
## restless 0.5589534 0.3479190 0.17381758 0.3391959 0.2932832
## sorestless 1.0000000 0.3547217 0.10398320 0.3999832 0.3596118
## depressed 0.3547217 1.0000000 0.13026042 0.3786149 0.3424609
## everythingeffort 0.1039832 0.1302604 1.00000000 0.0810859 0.1387534
## nothingcheerup 0.3999832 0.3786149 0.08108590 1.0000000 0.3980457
## worthless 0.3596118 0.3424609 0.13875344 0.3980457 1.0000000We find no strong correlations (>0.7) between an individual’s responses to the Kessler-10 questions. Thus, we cannot reliably use correlated variables to impute missing values. We must turn to another strategy. If individuals are unlikely to respond to the mental health questions in vastly different ways then we can say that the individual’s most frequently occuring response is the most likely value for the missing data. This implies a mode imputation where an individual’s most frequently occuring rating on the Kessler-10 scale will be used to impute missing values. We will calculate the mode in both waves since it is reasonable to believe that mental health can change over periods of time.
# write a function to calculate mode
# reference: https://stackoverflow.com/questions/2547402/how-to-find-the-statistical-mode
get_mode <- function(x) {
tab <- table(x)
tab <- tab[!is.na(names(tab))] # Remove NA values from table
if (length(tab) == 0) return(0) # Return 0 if all values in the vector are NA
return(names(tab)[which.max(tab)])
}
# Kessler-10 scale questions (mental health indicators)
kessler_columns <- c('tired', 'nervous', 'sonervous','hopeless', 'restless', 'sorestless',
'depressed', 'everythingeffort', 'nothingcheerup', 'worthless')
# impute the mode for each household member across waves
depression_imputed <- depression %>%
group_by(hhid, hhmid, wave) %>%
mutate(across(all_of(kessler_columns), ~ifelse(is.na(.), get_mode(na.omit(.)), .), .names = "imputed_{.col}")) %>%
ungroup()
# check that the imputation worked
colSums(is.na(depression_imputed))## wave hhid hhmid
## 0 0 0
## tired nervous sonervous
## 31 34 40
## hopeless restless sorestless
## 34 88 48
## depressed everythingeffort nothingcheerup
## 45 60 94
## worthless imputed_tired imputed_nervous
## 88 0 0
## imputed_sonervous imputed_hopeless imputed_restless
## 0 0 0
## imputed_sorestless imputed_depressed imputed_everythingeffort
## 0 0 0
## imputed_nothingcheerup imputed_worthless
## 0 0# We have now created new columns with the prefix "imputed_" to show that they
#are imputed and do not contain missing values
# calculate the Kessler score: sum of the imputed responses to the Kessler-10 questions
depression_imputed <- depression_imputed %>%
mutate(kessler_score = rowSums(dplyr::select(depression_imputed, starts_with("imputed_")), na.rm = TRUE)) %>%
dplyr::select(-tired, -nervous, -sonervous, -hopeless, -restless, -sorestless, -depressed, -everythingeffort, -nothingcheerup, -worthless)
# assign Kessler categories to the scores
depression_imputed <- depression_imputed %>%
mutate(kessler_category = case_when(
kessler_score <= 19 ~ "No Significant Depression",
kessler_score >= 20 & kessler_score <= 24 ~ "Mild Depression",
kessler_score >= 25 & kessler_score <= 29 ~ "Moderate Depression",
kessler_score >= 30 ~ "Severe Depression"
))Constructing a single dataset
Q7
At this point you have created three datasets: demographics, assets, and mental health. Please combine all three of these datasets to create a single dataset that you will use for data exploration and analysis. The unit of observation in this dataset should be an individual in a given survey round. (There should be at most two observations per individual, one for Wave 1 and another for Wave 2).
# first combine the demographics data with depression data by matching on wave,
#hhid, hhmid
temp_data <- inner_join(demographics_proxy, depression_imputed, by= c('wave','hhid', 'hhmid'))
# rename columns to match the asset data and remove the "imputed_" prefix
new_colnames <- c("Wave ID", "Household ID", "Household Member ID",
"tired", "nervous","sonervous","hopeless", "restless", "sorestless",
"depressed", "everythingeffort", "nothingcheerup", "worthless")
old_colnames <- c("wave" , "hhid", "hhmid", "imputed_tired", "imputed_nervous", "imputed_sonervous", "imputed_hopeless", "imputed_restless","imputed_sorestless", "imputed_depressed", "imputed_everythingeffort", "imputed_nothingcheerup", "imputed_worthless" )
temp_data <- temp_data %>% rename_with(~new_colnames, all_of(old_colnames))
# join the asset data by matching on wave and household id
panel_dta <- temp_data %>% inner_join(assets_hhwave, by= c("Wave ID", "Household ID"))
panel_csv <- write_csv(panel_dta, "panel_data.csv")Part II: Exploratory Analysis
Using Wave 1 data, conduct exploratory analysis to understand the relationship between depression and household and demographic characteristics among individuals in Ghana. Specifically, do the following:
Q1: Explore the relationship between depression and:
1.Household wealth, proxied by total asset value.
wave1_panel <- panel_dta %>% filter(`Wave ID` == 1) # filter data to only keep wave 1 obs
# visualize distribution of total asset value
eda_hist <- ggplot(wave1_panel, aes(x = `Total asset value`)) +
geom_histogram(bins = 50, fill = "blue", alpha = 0.6) +
labs(title = "Distribution of Total Asset Value",
x = "Total Asset Value",
y = "Count") +
theme_minimal()
ggsave("Distribution of Total Asset Value.jpg", plot = eda_hist, width = 8, height = 6, dpi = 300)
# scatter plot of depression and total asset value
eda_scatter <- ggplot(wave1_panel, aes(x = `Total asset value` , y = kessler_score)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", color = "red") +
scale_x_log10() + # Log scale to handle skewed asset values
labs(title = "Kessler Score vs. Total Asset Value",
x = "Total Asset Value",
y = "Kessler Score") +
theme_minimal()
ggsave("Scatter Plot - Depression and Marital Status.jpg", plot = eda_scatter, width = 8, height = 6, dpi = 300)## `geom_smooth()` using formula = 'y ~ x'eda_scatter## `geom_smooth()` using formula = 'y ~ x'
The total asset values in this dataset are skewed right with majority of
the values falling between 0 to 5000 GHS. There are some very high
values in the 10,000 to 15,000 GHS range indicating outliers. This
variable will be hard to work without a log transformation. Thus, we
will use the log values of total asset to explore the relationship
between household wealth and depression.
There is no clear linear relationship between asset value and depression. Majority values occur in the range of 500 to 5000 GHS in asset value and a Kessler score of 10 - 35. Low values of the Kessler score (<20), indicating mild to no symptoms of depression, occur more frequently in the 500 to 1000 GHS asset range while higher values of the Kessler scores (>30) are more concentrated in the 1000 to 1500 GHS asset value range. There is a smaller cluster of values in 75 to 500 GHS range with higher number of Kessler scores >20.
- A household or demographic characteristic that seems interesting to you.
I chose to explore the relationship between depression and marital status. Differences in marital status can produce vast differences in social status which can directly impact mental health. Married couples may enjoy better mental health due to social acceptance, companionship, and sharing of responsibility while separated or divorced individuals might face social exclusion and stigma resulting in poorer mental health. This effect can be even more pronounced by the gender.Since the treatment eventually invites heads of household and their spouses to attend group therapy, I am curious to know if there is a correlational relationship between marital status and depression. In this, I will also add gender to visualize the gendered impact on depression by marital status
unique(wave1_panel$maritalstatus)## <labelled<double>[8]>: 22 What is Name’s marital status?
## [1] 1 6 2 4 5 3 7 NA
##
## Labels:
## value label
## 1 Married
## 2 Consensual union
## 3 Separated
## 4 Divorced
## 5 Widowed
## 6 Never married
## 7 Betrothedlabels <- c("Married","Consensual union","Separated","Divorced","Widowed",
"Never married","Betrothed")
# box plot of depression and marital status
# reference: https://r-charts.com/ggplot2/axis/
eda_boxplot <- ggplot(wave1_panel, aes(x = as.factor(maritalstatus), y = kessler_score)) +
geom_boxplot(fill = "dodgerblue1",
colour = "black",
alpha = 0.5,
outlier.colour = "tomato2") +
labs(title = "Distribution of Kessler Score by Marital Status", x = "Marital Status", y = "Kessler Score") +
theme_minimal() +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_x_discrete(label = labels)
ggsave("Box plot - Depression and Marital Status.jpg", plot = eda_boxplot, width = 8, height = 6, dpi = 300)
eda_boxplot
Evaluating the RCT
Using Wave 2 data to measure outcomes, answer the following questions, explaining any decisions and assumptions you make, and interpret your results. There is no need for you to address the validity of the random assignment of the intervention.
Q2:
Were the GT sessions effective at reducing depression?
\(\color{red}{\text{SOLUTION:}}\)
To answer this question, we first need to pick an outcome measuring depression. There are several depression variables (“depressed”, “Kessler score”, “Kessler category”, etc) that we could consider as the primary outcome variable for measuring changes between the treated and the untreated. In order to have a continuous dependent variable in the regression I chose the Kessler score as the outcome variable measuring depression. Since we are not concerned with the validity of the random assignment, we assume that the households that got the treatment were identical to the households that did not get the treatment in all respects except their treatment condition. Thus, after controlling for observable characteristics like age, gender, marital status, household size, etc., we expect that any changes in the Kessler score would be caused by the treatment i.e. receiving group therapy.
wave2_panel <- panel_dta %>% filter(`Wave ID` == 2)
wave2_paneldepression_reg1 <-
lm(kessler_score ~ treat_hh + gender + age + maritalstatus + relationship + hhsize, data = wave2_panel)
stargazer(depression_reg1,
type = "text",
header = FALSE,
title = "Impact of Group Therapy on Depression",
column.labels = c("$depression_reg1.$"),
colnames = FALSE,
model.numbers = FALSE,
df = FALSE)##
## Impact of Group Therapy on Depression
## ===============================================
## Dependent variable:
## ---------------------------
## kessler_score
## depression
## -----------------------------------------------
## treat_hh -0.515***
## (0.135)
##
## gender 0.175***
## (0.049)
##
## age 0.032***
## (0.003)
##
## maritalstatus 0.229***
## (0.051)
##
## relationship -0.557**
## (0.225)
##
## hhsize 0.172***
## (0.028)
##
## Constant 14.945***
## (0.366)
##
## -----------------------------------------------
## Observations 6,382
## R2 0.033
## Adjusted R2 0.032
## Residual Std. Error 5.384
## F Statistic 36.605***
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01# Plot treatment effects
ggplot(wave2_panel, aes(x = as.factor(treat_hh), y = kessler_score)) +
geom_boxplot() +
labs(title = "Effect of Treatment on Kessler Score", x = "Treatment Group", y = "Kessler Score") +
theme_minimal()
Q3
Did the effect of GT sessions on depression differ by gender? Perform a linear regression of the Kessler Score on:
Woman (binary variable)
Treated Household (binary variable)
Interaction term: Treated Household
unique(wave2_panel$gender) # check which values gender takes ## <labelled<double>[2]>: Gender
## [1] 1 5
##
## Labels:
## value label
## 1 Male
## 5 Female# create a binary variable named woman
wave2_panel <- wave2_panel %>% mutate(woman = ifelse(gender == 5, 1, 0))depression_reg2 <-
lm(kessler_score ~ treat_hh + woman + treat_hh * woman, data = wave2_panel)
stargazer(depression_reg2,
type = "text",
header = FALSE,
title = "Impact of Group Therapy on Depression by Gender",
column.labels = c("$depression_reg2.$"),
colnames = FALSE,
model.numbers = FALSE,
df = FALSE)##
## Impact of Group Therapy on Depression by Gender
## ===============================================
## Dependent variable:
## ---------------------------
## kessler_score
## depression
## -----------------------------------------------
## treat_hh -0.132
## (0.207)
##
## woman 0.796***
## (0.195)
##
## treat_hh:woman -0.706**
## (0.276)
##
## Constant 17.063***
## (0.147)
##
## -----------------------------------------------
## Observations 6,382
## R2 0.005
## Adjusted R2 0.005
## Residual Std. Error 5.461
## F Statistic 10.648***
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01