Homework 2
2025-05-14
1 Introduction
Depression is a global public health priority. According to the World Health Organization and the World Bank depression accounts for a full 10% of the total nonfatal burden of disease worldwide (World Health Organization & World Bank, 2016). Understanding the factors that contribute to depression is critical for developing effective interventions and policies to mitigate its impact on individuals and society. This study aims to explore the relationship between depression and key variables, including education level, gender, frequency of vegetable consumption and physical activity. These variables were selected as they encompass both socio-demographic and lifestyle influences known to affect mental health. Focusing on Ireland, a country that has participated in all rounds of the European Social Survey (ESS), this study draws on data from over 2,000 respondents to provide insights into the factors associated with depression.
Research Question:
How do education level, gender, frequency of vegetable consumption and physical activity, influence the prevalence or severity of depression among the population of Iceland?
Dependent Variable:
• Depression
Independent Variables:
• Education level
• Gender
• Frequency of vegetable consumption
• Physical activity
As a first step we load all libraries needed and load the data we wil be working with. The ESS11.sav data from the European Social Survey.
library(foreign)
library(ltm)## Loading required package: MASS## Loading required package: msm## Loading required package: polycorsetwd("/Users/lenapape/Desktop/Uni /Quantative Research R")
df = read.spss("ESS11.sav", to.data.frame = T)
#names(df)
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)In the following section we created the depression scale using the
eight factors of the CES-D8 Depression Scale, which assesses depressive
symptoms based on responses to 8 items:
1. ‘fltdpr’: how often have you felt depressed in the past week
2. ‘flteeff’: how often did everything you did in the past week feel as
an effort?
3. ‘slprl’: Sleep was restless, how often past week
4. ‘wrhpp’: Were happy, how often past week
5. ‘fltlnl’: Felt lonely, how often past week
6. ‘enjlf’: Enjoyed life, how often past week
7. ‘fltsd’: Felt sad, how often past week
8. ‘cldgng’: Could not get going, how often past week
All responses are converted into a numeric scala with values ranging from 1 to 4.
#convert to numbers 1-4
df$d20 = as.numeric(df$fltdpr)
df$d21 = as.numeric(df$flteeff)
df$d22 = as.numeric(df$slprl)
df$d23 = as.numeric(df$wrhpp)
df$d24 = as.numeric(df$fltlnl)
df$d25 = as.numeric(df$enjlf)
df$d26 = as.numeric(df$fltsd)
df$d27 = as.numeric(df$cldgng)
df$d23 = 5 - df$d23
df$d25 = 5 - df$d25
cronbach.alpha(df[,c("d20", "d21", "d22", "d23", "d24", "d25", "d26", "d27")], na.rm=T)##
## Cronbach's alpha for the 'df[, c("d20", "d21", "d22", "d23", "d24", "d25", "d26", "d27")]' data-set
##
## Items: 8
## Sample units: 40156
## alpha: 0.823The following table presents the distribution of answers among the eight different variables in our depression scale. The chart below visualizes these numbers further.
library(likert)
library(kableExtra)
library(knitr)
vnames = c("fltdpr", "flteeff", "slprl", "wrhpp", "fltlnl", "enjlf", "fltsd", "cldgng")
likert_df = df[,vnames]
likert_table = likert(likert_df)$results
likert_numeric_df = as.data.frame(lapply((df[,vnames]), as.numeric))
likert_table$Mean = unlist(lapply((likert_numeric_df[,vnames]), mean, na.rm=T)) # ... and append new columns to the data frame
likert_table$Count = unlist(lapply((likert_numeric_df[,vnames]), function (x) sum(!is.na(x))))
likert_table$Item = c(
d20="How much of the time during the past week did you feel depressed?",
d21="How much of the time during the past week did you feel that everything you did was an effort?",
d22="How much of the time during the past week was your sleep restless",
d23="How much of the time during the past week you were happy?",
d24="How much of the time during the past week did you feel lonely?",
d25="How much of the time during the past week did you enjoyed life?",
d26="How much of the time during the past week did you feel sad?",
d27="How much of the time during the past week did you feel you could not get going?"
)
# round all percentage values to 1 decimal digit
likert_table[,2:6] = round(likert_table[,2:6],1) # round columns 2 to 6 to 1 digit
# round means to 3 decimal digits
likert_table[,7] = round(likert_table[,7],3) # round column number 7 to 3 digits
kable_styling(kable(likert_table,
caption = "Distribution of answers regarding depression (ESS round 11, all countries)"
)
)| Item | None or almost none of the time | Some of the time | Most of the time | All or almost all of the time | Mean | Count |
|---|---|---|---|---|---|---|
| How much of the time during the past week did you feel depressed? | 64.9 | 29.1 | 4.6 | 1.5 | 1.4 | 39981 |
| How much of the time during the past week did you feel that everything you did was an effort? | 48.4 | 38.4 | 9.8 | 3.4 | 1.7 | 39983 |
| How much of the time during the past week was your sleep restless | 43.9 | 39.9 | 11.6 | 4.6 | 1.8 | 40017 |
| How much of the time during the past week you were happy? | 4.0 | 23.5 | 48.9 | 23.6 | 2.9 | 39890 |
| How much of the time during the past week did you feel lonely? | 68.1 | 24.3 | 5.3 | 2.3 | 1.4 | 39983 |
| How much of the time during the past week did you enjoyed life? | 5.3 | 24.8 | 44.8 | 25.0 | 2.9 | 39878 |
| How much of the time during the past week did you feel sad? | 52.5 | 41.1 | 4.9 | 1.6 | 1.6 | 39981 |
| How much of the time during the past week did you feel you could not get going? | 55.7 | 36.1 | 6.2 | 2.0 | 1.5 | 39949 |
# create basic plot (code also valid)
plot(likert(summary=likert_table[,1:6]))
#install.packages("psych")
library(psych)
alpha_result = alpha(df[, c("d20", "d21", "d22", "d23", "d24", "d25", "d26", "d27")])
alpha_rnd = round(alpha_result$total[["raw_alpha"]], 3)
# alpha_rndTo show if the scale demonstartes high internal reliability we calculated the cronbachs alpha value. In our case we received an alpha value of 0.824, this shows a high reliability because it is above 0.7.
To compute an overall depression score, we averaged the values of the
8 items (d20–d27) for each respondent. This new variable
depression reflects the mean score across the items, and we
used summary() to examine its distribution.
df$depression = rowSums(df[,c("d20", "d21", "d22", "d23", "d24", "d25", "d26", "d27")]) / 8
summary(df$depression)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 1.375 1.625 1.695 2.000 4.000 799Furthermore, to show the distribution of the results of the depression scala in a histogram:
library(kableExtra)
library(knitr)
hist(df$depression,
breaks=8,
col = "lightblue",
border = "white",
main = paste("Histogramm of Depression Score"),
xlab = "Depression Score")
table_depression = as.data.frame(table(df$depression))
#kable( table_depression,
# col.names = c("Varianz","Frequency"),
#caption = "Distribution of Depression")
# scroll_box(
# kable_styling(
# kable(table_depression,
# col.names = c("Varianz","Frequency"),
# caption = "Distribution of Depression"
# ),
# font_size = 13, bootstrap_options = c("hover", "condensed")
# ),
# height="200px")
In the next section we visualize which countries gathered data on the
depression scale. And with these results we selected the country of
Ireland.
# Subset the data to a single country; Ireland
# df$cntry
library(kableExtra)
library(knitr)
table_country = table(droplevels(df$cntry))
table_country = as.data.frame(table_country)
scroll_box(
kable_styling(
kable(table_country,
col.names = c("Country","Data"),
caption = "Data Entries per Country"
),
full_width = F, font_size = 13, bootstrap_options = c("hover", "condensed")
),
height="300px")| Country | Data |
|---|---|
| Austria | 2354 |
| Belgium | 1594 |
| Switzerland | 1384 |
| Cyprus | 685 |
| Germany | 2420 |
| Spain | 1844 |
| Finland | 1563 |
| France | 1771 |
| United Kingdom | 1684 |
| Greece | 2757 |
| Croatia | 1563 |
| Hungary | 2118 |
| Ireland | 2017 |
| Iceland | 842 |
| Italy | 2865 |
| Lithuania | 1365 |
| Netherlands | 1695 |
| Norway | 1337 |
| Poland | 1442 |
| Portugal | 1373 |
| Serbia | 1563 |
| Sweden | 1230 |
| Slovenia | 1248 |
| Slovakia | 1442 |
Choose Country of Ireland and create new data frame only consisting of data collected in Ireland
df_irl= df[df$cntry == "Ireland",]
#df_irl
Since we will only be analysing the depression scale of Ireland we
repeat the computation of the frequency distribution, but only
exclusivley for data from Ireland.
table_irl = as.data.frame(table(df_irl$depression))
scroll_box(
kable_styling(
kable(table_irl,
col.names = c("Varianz","Frequency"),
caption = "Distribution of Depression in Ireland"
),
full_width = F, font_size = 13, bootstrap_options = c("hover", "condensed")
),
height="200px")| Varianz | Frequency |
|---|---|
| 1 | 269 |
| 1.125 | 207 |
| 1.25 | 236 |
| 1.375 | 210 |
| 1.5 | 203 |
| 1.625 | 159 |
| 1.75 | 139 |
| 1.875 | 125 |
| 2 | 113 |
| 2.125 | 79 |
| 2.25 | 88 |
| 2.375 | 45 |
| 2.5 | 19 |
| 2.625 | 19 |
| 2.75 | 12 |
| 2.875 | 8 |
| 3 | 10 |
| 3.125 | 4 |
| 3.25 | 5 |
| 3.375 | 2 |
| 3.5 | 2 |
| 3.875 | 1 |
hist(df_irl$depression, breaks = seq(1, 4, by = 0.25),
angle = 45,
col = "lightblue",
border = "white",
main = paste("Histogramm of Depression Score in Ireland"),
xlab = "Depression Score")
2 Statistical Description
For a brief overview we calculated how many of our participants in the data set were female and male as well as the age distribution. ### Gender Distribution{.unnumbered}
table(df_irl$gndr)##
## Male Female
## 906 1111sum(df_irl$gndr == 'Male')## [1] 906sum(df_irl$gndr == 'Female')## [1] 1111A total of 906 males and 1111 females are in the dataset. And below we computed a histogram showing the age distribution among the survey participants.
Age Distribution
#head(df_irl$agea)
#summary(df_irl$agea)
df_irl$agea = as.numeric(as.character(df_irl$agea))
hist(df_irl$agea,
main = "Age Distribution",
xlab = "Age",
col = "lightblue",
border = "white")
3 Hypothesis
3.1 Hypothesis 1 - Education Level:
The relationship between education levels and the prevalence of depression is well-documented, indicating significant disparities in mental health outcomes. Research consistently shows that individuals with lower educational levels have a 19,1% higher likelihood to experience major depressive disorder (MDD) between the age of 18 and 65 compared to those with higher education (Lepe et al., 2022). Moreover, a meta-analysis from 2022 shows a significant prevalence of depressive symptoms, like anxiety, among adolescents with low education and income (Kempfer et al., 2022).
H1: Hypothesis one states that a lower socioeconomic status (education level) is associated with higher levels of depression. We assume that lower education levels will be positively associated with higher CES-D8 depression scores.
# BIVARIATE ANALYSIS AND OPERATIONALIZATION
# Hypothesis 1: Socioeconomic Status: Lower socioeconomic status (education level) is associated
# with higher levels of depression
# recode education levels into 3 groups
df_irl$edu = factor(NA, levels = c("low", "medium", "high")) # variables created as factors with levels (low, medium, high)
# look up original values
# table(df_irl$eisced)
Next, we recoded education into three groups, aiming to distribute the
data as evenly as possible across the categories.
kable_styling(
kable(as.data.frame(table(df_irl$edu)),
col.names = c("Level of Education","Distribution"),
caption = "Distribution Education Level in Ireland"
),
full_width = F, font_size = 13, bootstrap_options = c("hover", "condensed")
)| Level of Education | Distribution |
|---|---|
| low | 668 |
| medium | 756 |
| high | 583 |
library(knitr)
library(kableExtra)
# Create the data frame
edu_table <- data.frame(
`Education Level` = c(
"Not possible to harmonise into ES-ISCED",
"ES-ISCED I , less than lower secondary",
"ES-ISCED II, lower secondary",
"ES-ISCED IIIb, lower tier upper secondary",
"ES-ISCED IIIa, upper tier upper secondary",
"ES-ISCED IV, advanced vocational, sub-degree",
"ES-ISCED V1, lower tertiary education, BA level",
"ES-ISCED V2, higher tertiary education, >= MA level",
"Other"
),
Low = c(0, 197, 376, 95, 0, 0, 0, 0, 0),
Medium = c(0, 0, 0, 0, 320, 436, 0, 0, 0),
High = c(0, 0, 0, 0, 0, 0, 298, 285, 0)
)
# Create the nicely formatted table
kable_styling(
kable(edu_table,
col.names = c("Education Level", "Low", "Medium", "High"),
caption = "Mapping of ES-ISCED Categories to Education Groups"),
full_width = FALSE,
font_size = 13,
bootstrap_options = c("hover", "condensed")
)| Education Level | Low | Medium | High |
|---|---|---|---|
| Not possible to harmonise into ES-ISCED | 0 | 0 | 0 |
| ES-ISCED I , less than lower secondary | 197 | 0 | 0 |
| ES-ISCED II, lower secondary | 376 | 0 | 0 |
| ES-ISCED IIIb, lower tier upper secondary | 95 | 0 | 0 |
| ES-ISCED IIIa, upper tier upper secondary | 0 | 320 | 0 |
| ES-ISCED IV, advanced vocational, sub-degree | 0 | 436 | 0 |
| ES-ISCED V1, lower tertiary education, BA level | 0 | 0 | 298 |
| ES-ISCED V2, higher tertiary education, >= MA level | 0 | 0 | 285 |
| Other | 0 | 0 | 0 |
prop.table(table(df_irl$edu))##
## low medium high
## 0.3328351 0.3766816 0.2904833# Example: Create weights based on edu distribution
edu_freq <- table(df_irl$edu)
edu_weights <- 1 / as.numeric(edu_freq)
names(edu_weights) <- names(edu_freq)
# Assign weights to each row based on their edu level
df_irl$weight <- edu_weights[as.character(df_irl$edu)]#check significance
kw_test = kruskal.test(depression ~ edu, data=df_irl)
kw_test##
## Kruskal-Wallis rank sum test
##
## data: depression by edu
## Kruskal-Wallis chi-squared = 31.065, df = 2, p-value = 1.796e-07p_edu = round(kw_test$p.value, 10)
The Kruskal-Wallis test for this variable yielded a p-value of
1.796^{-7}, indicating a statistically significant difference between
groups (p < 0.05).
by(df_irl$depression, df_irl$edu, mean, na.rm=T)## df_irl$edu: low
## [1] 1.638351
## ------------------------------------------------------------
## df_irl$edu: medium
## [1] 1.567425
## ------------------------------------------------------------
## df_irl$edu: high
## [1] 1.468805model_edu = lm(depression ~ edu, data = df_irl)
summary(model_edu)##
## Call:
## lm(formula = depression ~ edu, data = df_irl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63835 -0.34381 -0.09381 0.30757 2.40619
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.63835 0.01817 90.179 < 2e-16 ***
## edumedium -0.07093 0.02490 -2.849 0.00443 **
## eduhigh -0.16955 0.02660 -6.374 2.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4618 on 1944 degrees of freedom
## (70 observations deleted due to missingness)
## Multiple R-squared: 0.02054, Adjusted R-squared: 0.01954
## F-statistic: 20.39 on 2 and 1944 DF, p-value: 1.726e-093.2 Hypothesis 2 - Gender:
Myrna Weissman pioneered the field of gender-based research in depression in the 1970s, highlighting a significant gender disparity. She observed that among clinical and community samples, the prevalence of depression was twice as high among women as it was among men (Weissman & Klerman, 1977). After this landmark a lot of other research was conducted in this field (e.g. Bebbington et al., 1998; Ferrari et al., 2012) However, the most observed ratio of female to male cases typically ranges from to one, indicating a clear disparity in the prevalence of the condition between the two genders.
H2: Women report higher levels of depression compared to men. → Women will score higher on the CES-D8 scale compared to men.
# variable: 'gndr'
#df_irl$gndr
# test significance
wk_test = wilcox.test(depression ~ gndr, data = df_irl)
wk_test##
## Wilcoxon rank sum test with continuity correction
##
## data: depression by gndr
## W = 443514, p-value = 0.01777
## alternative hypothesis: true location shift is not equal to 0p_gndr = round(wk_test$p.value, 4)
by(df_irl$depression, df_irl$gndr, mean, na.rm=T)## df_irl$gndr: Male
## [1] 1.540575
## ------------------------------------------------------------
## df_irl$gndr: Female
## [1] 1.580316# create female dummy
table(df_irl$gndr)##
## Male Female
## 906 1111df_irl$female = NA #intiialize variable with NA
df_irl$female[df_irl$gndr=="Male"] = 0 # if variable is female then gender == 0
df_irl$female[df_irl$gndr=="Female"] = 1 # if variable is male then gender == 1
# check
table(df_irl$female)##
## 0 1
## 906 1111table(df_irl$gndr, df_irl$female)##
## 0 1
## Male 906 0
## Female 0 1111lm(depression ~ female, data = df_irl)##
## Call:
## lm(formula = depression ~ female, data = df_irl)
##
## Coefficients:
## (Intercept) female
## 1.54058 0.03974lm(depression ~ gndr, data = df_irl)##
## Call:
## lm(formula = depression ~ gndr, data = df_irl)
##
## Coefficients:
## (Intercept) gndrFemale
## 1.54058 0.03974
The Wilcox Test showed statisticak sugnificance becaus the p value lied
at 0.0178, which is below 0.05. The linear regression model also
indicates that depression is more prevelant among females rather than
males.
3.3 Hypothesis 3 - Vegetable Consumption:
The relationship between vegetable consumption and depression is complex, with a variety of studies in this area. A national survey among the Canadian population found out that people with a high intake of daily vegetables and fruits had a lower risk of major depressive disorder (MDD) (McMartin et al., 2013). Moreover, a systematic review indicated that a higher intake of fruits and vegetables has a positive effect on women’s mental health (Guzek et al., 2022). These findings align with the recommendations for the prevention of MDD. This study identified five dietary recommendations for prevention including increased consumption of fruits, vegetables, legumes, wholegrain cereals, nuts and seeds (Opie et al., 2017).
H3: Individuals who consume vegetables less frequently have higher levels of depression compared to those who consume them regularly. → Lower frequency of vegetable consumption will correspond to higher CES-D8 scores
# variable = eatveg
#df_irl$eatveg
levels(df_irl[,"eatveg"])## [1] "Three times or more a day"
## [2] "Twice a day"
## [3] "Once a day"
## [4] "Less than once a day but at least 4 times a week"
## [5] "Less than 4 times a week but at least once a week"
## [6] "Less than once a week"
## [7] "Never"
As with the education levels we recoded the vegetable intake levels into
three groups, aiming to distribute the data as evenly as possible across
the categories.
# change level names
df_irl$veggis = factor(NA, levels = c("multiple daily", "daily", "irregular"))
df_irl$veggis[df_irl$eatveg == "Three times or more a day"] = "multiple daily"
df_irl$veggis[df_irl$eatveg == "Twice a day"] = "multiple daily"
df_irl$veggis[df_irl$eatveg == "Once a day"] = "daily"
df_irl$veggis[df_irl$eatveg == "Less than once a day but at least 4 times a week"] = "irregular"
df_irl$veggis[df_irl$eatveg == "Less than 4 times a week but at least once a week"] = "irregular"
df_irl$veggis[df_irl$eatveg == "Less than once a week"] = "irregular"
df_irl$veggis[df_irl$eatveg == "Never"] = "irregular"
veg_table = as.data.frame(table(df_irl$veggis))
# veg_table
library(kableExtra)
library(knitr)
kable_styling(
kable(veg_table,
col.names = c("Vegetable Intake", "Frequency"),
caption = "3 levels of vegetable intake"
),
full_width = F, font_size = 13, bootstrap_options = c("hover", "condensed")
)| Vegetable Intake | Frequency |
|---|---|
| multiple daily | 557 |
| daily | 1052 |
| irregular | 407 |
#check significance
kt_veg = kruskal.test(depression ~ veggis, data=df_irl)
p_veg = kt_veg$p.value
by(df_irl$depression, df_irl$veggis, mean, na.rm=T)## df_irl$veggis: multiple daily
## [1] 1.561806
## ------------------------------------------------------------
## df_irl$veggis: daily
## [1] 1.527275
## ------------------------------------------------------------
## df_irl$veggis: irregular
## [1] 1.654898The Kruskal-Wallis Test resulted in a p value of 1.39^{-4} showing statistical significance. But the comparison of all 3 levels didin’t show irregular vegetable consumption having a clear effect. Since vegetable consumption of multiple daily had a higher value than daily consumption.
3.4 Hypothesis 4 - Physical Activity
Engagement in physical activity (PA) is significantly associated with lower levels of depression, according to several studies. A recently published cross-sectional study involving 58,455 adults found a clear link between lower levels of PA and increased depressive symptoms (Queiroga et al., 2025). Moreover, a multitude of studies have demonstrated the varying levels of PA can effectively diminish the likelihood of developing depression (e.g. Schuch et al., 2018). This is also reported by the National Health and Nutrition Examination Survey from 2007 to 2018, this cross-sectional analysis indicated that participants without depressive symptoms engage in significantly more PA (10.87 hours/week) than those with depressive symptoms (8.82 hours/week) (Bopari et al., 2024).
H4: Individuals who engage in physical activity on fewer days in the last week report higher levels of depression compared to those who are more physically active. → Fewer days of physical activity will correlate with higher CES-D8 scores.
We conducted a one-way ANOVA to examine whether the number of days an
individual engaged in sports or physical activity in the past 7 days
(dosprt) is associated with differences in depression
levels. Using the dataset df_irl, the model tested for mean
differences in depression scores across the 8 levels of
dosprt (ranging from 0 to 7 days).
# variable = dosprt : "Do sports or other physical activity, how many of last 7 days"
#df_irl$dosprt
# test significance
anova_m = aov(depression ~ dosprt, data = df_irl)
summary(anova_m)## Df Sum Sq Mean Sq F value Pr(>F)
## dosprt 7 22.9 3.276 15.9 <2e-16 ***
## Residuals 1945 400.6 0.206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 64 observations deleted due to missingnessp_value_sprts = summary(anova_m)[[1]][["Pr(>F)"]][1]The results indicated a statistically significant effect of physical activity on depression, F(7, 1945) = 15.9, p = 1.9649655^{-20}. This suggests that the number of physically active days is significantly related to reported depression levels. Note that 64 observations were excluded due to missing values.
# make dosprt variable numeric
df_irl$dosprt_n = as.numeric(df_irl$dosprt) -1
# pairwise correlation
# cor(df_irl[,c("depression", "dosprt_n")], use="complete.obs")
# results show negative value --> "the more sports, the less depressed"
library(knitr)
library(kableExtra)
# Calculate correlation matrix
cor_matrix <- cor(df_irl[, c("depression", "dosprt_n")], use = "complete.obs")
# Format as a data frame for table display
cor_df <- as.data.frame(round(cor_matrix, 2))
# Optional: move row names into a column
#cor_df$Variable <- rownames(cor_df)
#cor_df <- cor_df[, c("Variable", "depression", "dosprt_n")]
# Create the nicely styled table
kable_styling(
kable(cor_df,
caption = "Correlation Between Depression and Sport Frequency"),
full_width = FALSE,
font_size = 13,
bootstrap_options = c("hover", "condensed")
)| depression | dosprt_n | |
|---|---|---|
| depression | 1.00 | -0.23 |
| dosprt_n | -0.23 | 1.00 |
# linear regression model
# lm(depression ~ dosprt, data = df_irl, na.rm=T)
# Required packages
library(broom)
library(knitr)
library(kableExtra)
library(dplyr)
# Run linear regression model
model_sport <- lm(depression ~ dosprt, data = df_irl)
# Tidy the model output (with confidence intervals)
tidy_model <- tidy(model_sport, conf.int = TRUE)
# Round results for readability
tidy_model <- tidy_model %>%
mutate(across(where(is.numeric), ~ round(., 3)))
# Create the nicely formatted table
kable_styling(
kable(tidy_model,
col.names = c("Term", "Estimate", "Std. Error", "t value", "p value", "CI Lower", "CI Upper"),
caption = "Linear Regression: Effect of Sport Frequency on Depression Score"),
full_width = FALSE,
font_size = 13,
bootstrap_options = c("hover", "condensed")
)| Term | Estimate | Std. Error | t value | p value | CI Lower | CI Upper |
|---|---|---|---|---|---|---|
| (Intercept) | 1.712 | 0.021 | 81.151 | 0.000 | 1.670 | 1.753 |
| dosprt1 | -0.074 | 0.050 | -1.490 | 0.136 | -0.172 | 0.023 |
| dosprt2 | -0.075 | 0.039 | -1.903 | 0.057 | -0.152 | 0.002 |
| dosprt3 | -0.168 | 0.035 | -4.774 | 0.000 | -0.237 | -0.099 |
| dosprt4 | -0.153 | 0.038 | -3.993 | 0.000 | -0.227 | -0.078 |
| dosprt5 | -0.227 | 0.037 | -6.123 | 0.000 | -0.300 | -0.154 |
| dosprt6 | -0.208 | 0.056 | -3.753 | 0.000 | -0.317 | -0.100 |
| dosprt7 | -0.295 | 0.030 | -9.751 | 0.000 | -0.354 | -0.236 |
The Test yielded the following results: more physical activity (doing
sports on 3 - 7 days in a week) is linked to lower depression but lower
levels of physical activity (doing sports on 1 - 2 days in a week) don’t
show a significant difference from doing no sports at all
4 Multivariate Modelling
For a final overview we present a linear regression model with all variables included. The graph visualizes the results and in the table we present the values of each independent variable. This first regression model is unweighted.
#######################
# MULTIVARIATE MODELLING
# When hypothesis buidling and bivariate analysis is finalized,
# put everything together in a multiple regression model
# Unweighted model
model_unweighted = lm(depression ~ edu + female + veggis + dosprt, data = df_irl)
# summary(model_unweighted)
plot(model_unweighted, which =1)
# Model unweighted in tidy table
# Tidy the model output
model_tidy_unweighted = tidy(model_unweighted) # 'model' is your lm() object
# Round for readability
model_tidy_unweighted = model_tidy_unweighted %>%
mutate(across(where(is.numeric), ~ round(., 4)))
# Present it
kable(model_tidy_unweighted, caption = "Regression Output Unweighted: Coefficients of Depression Model")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 1.7440 | 0.0335 | 52.0206 | 0.0000 |
| edumedium | -0.0429 | 0.0250 | -1.7139 | 0.0867 |
| eduhigh | -0.1221 | 0.0273 | -4.4658 | 0.0000 |
| female | 0.0514 | 0.0208 | 2.4732 | 0.0135 |
| veggisdaily | -0.0679 | 0.0244 | -2.7790 | 0.0055 |
| veggisirregular | 0.0290 | 0.0311 | 0.9322 | 0.3513 |
| dosprt1 | -0.0452 | 0.0501 | -0.9029 | 0.3667 |
| dosprt2 | -0.0611 | 0.0395 | -1.5487 | 0.1216 |
| dosprt3 | -0.1440 | 0.0357 | -4.0347 | 0.0001 |
| dosprt4 | -0.1292 | 0.0387 | -3.3397 | 0.0009 |
| dosprt5 | -0.1911 | 0.0377 | -5.0639 | 0.0000 |
| dosprt6 | -0.1748 | 0.0559 | -3.1263 | 0.0018 |
| dosprt7 | -0.2641 | 0.0312 | -8.4768 | 0.0000 |
Next we present the weighted model using the post-stratification weight which was predefined by the European Social Survey (ESS) for their survey data. Since we are only looking at data from one country, Ireland in this case, we only require the pspwght weight. Below the results are presented in a graph and a table.
# Weighted model
model_weighted <- lm(depression ~ edu + female + veggis + dosprt, data = df_irl, weights = pspwght)
plot(model_weighted, which =1)
# summary(model_weighted)
# Install if needed
# install.packages("broom")
library(broom)
library(knitr)
# Model weighted in tidy table
# Tidy the model output
model_tidy_weighted <- tidy(model_weighted) # 'model' is your lm() object
# Round for readability
model_tidy_weighted <- model_tidy_weighted %>%
mutate(across(where(is.numeric), ~ round(., 4)))
# Present it
kable(model_tidy_weighted, caption = "Regression Output Weighted: Coefficients of Depression Model", full_width = F)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 1.7068 | 0.0350 | 48.8163 | 0.0000 |
| edumedium | -0.0128 | 0.0252 | -0.5084 | 0.6113 |
| eduhigh | -0.0780 | 0.0272 | -2.8616 | 0.0043 |
| female | 0.0396 | 0.0203 | 1.9535 | 0.0509 |
| veggisdaily | -0.0896 | 0.0236 | -3.7949 | 0.0002 |
| veggisirregular | -0.0302 | 0.0302 | -1.0005 | 0.3172 |
| dosprt1 | -0.0663 | 0.0493 | -1.3444 | 0.1790 |
| dosprt2 | -0.0538 | 0.0404 | -1.3312 | 0.1833 |
| dosprt3 | -0.0780 | 0.0355 | -2.1991 | 0.0280 |
| dosprt4 | -0.0756 | 0.0385 | -1.9623 | 0.0499 |
| dosprt5 | -0.1576 | 0.0372 | -4.2354 | 0.0000 |
| dosprt6 | -0.1657 | 0.0523 | -3.1709 | 0.0015 |
| dosprt7 | -0.2346 | 0.0322 | -7.2861 | 0.0000 |
The comparison only show little differences in values.
The graph below presents the estimated effects of all independent
variables on the depression score, including confidence intervals for
each regression coefficient, with the weighted model.
# visualize the model
#install.packages("ggplot2")
#install.packages("broom")
library(ggplot2)
library(broom)
tidy_model = tidy(model_weighted, conf.int = TRUE)
ggplot(tidy_model, aes(x = estimate, y = term)) +
geom_point() +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2) +
labs(title = "Regression Coefficients for Depression Model", x = "Coefficient Estimate", y = "Predictor") +
theme_minimal()