Social Anxiety
2025-03-19
Yong-Quan Su
The dataset contains 11000 individuals from different demographic background reporting their lifestyle, health & mental status, mental health history, life events, and their anxiety level. The analysis aims to find out whether there is any difference in caffeine intake at each reported anxiety level.
Import data
df <- read.table("./enhanced_anxiety_dataset.csv", header = TRUE, sep = ",", dec = ".")
head(df)## Age Gender Occupation Sleep.Hours Physical.Activity..hrs.week.
## 1 29 Female Artist 6.0 2.7
## 2 46 Other Nurse 6.2 5.7
## 3 64 Male Other 5.0 3.7
## 4 20 Female Scientist 5.8 2.8
## 5 49 Female Other 8.2 2.3
## 6 53 Male Other 6.4 6.5
## Caffeine.Intake..mg.day. Alcohol.Consumption..drinks.week. Smoking
## 1 181 10 Yes
## 2 200 8 Yes
## 3 117 4 No
## 4 360 6 Yes
## 5 247 4 Yes
## 6 235 2 No
## Family.History.of.Anxiety Stress.Level..1.10. Heart.Rate..bpm.
## 1 No 10 114
## 2 Yes 1 62
## 3 Yes 1 91
## 4 No 4 86
## 5 No 1 98
## 6 No 9 84
## Breathing.Rate..breaths.min. Sweating.Level..1.5. Dizziness Medication
## 1 14 4 No Yes
## 2 23 2 Yes No
## 3 28 3 No No
## 4 17 3 No No
## 5 19 4 Yes Yes
## 6 14 3 No Yes
## Therapy.Sessions..per.month. Recent.Major.Life.Event Diet.Quality..1.10.
## 1 3 Yes 7
## 2 2 No 8
## 3 1 Yes 1
## 4 0 No 1
## 5 1 No 3
## 6 2 Yes 5
## Anxiety.Level..1.10.
## 1 5
## 2 3
## 3 1
## 4 2
## 5 1
## 6 4Data description
Unit of observation: one individual
Sample size: 11000 observations
Variables and its units of measurement:
- Age: the age of each observation (numeric: 1 year)
- Gender: the gender of each observation (factorial: male/female/other)
- Occupation: the occupation of each observation (factorial: different occupation groups)
- Sleep.Hours: average sleep hours of each observation (numeric: 1 hour)
- Physical.Activity..hrs.week.: average hours each observation committing to physical activity (numeric: 1 hour)
- Caffeine.Intake..mg.day.: average caffeine intake of each observation per day (numeric: 1 mg)
- Alcohol.Consumption..drinks.week.: average alcohol consumption of each observation per week (numeric: 1 unit of alcohol)
- Smoking: smoking habit of each observation (factorial: yes/no)
- Family.History.of.Anxiety: history of anxiety within each observation’s family (factorial: yes/no)
- Stress.Level..1.10.: self-reported stress level of each observation (factorial: 1-10)
- Heart.Rate..bpm.: heart rate of each observation (numeric: 1 beat per minute)
- Breathing.Rate..breaths.min.: breathing rate of
each observation per minute (numeric: 1 breath per minute)
- Sweating.Level..1.5.: self-reported sweating level of each observation (factorial: 1-5)
- Dizziness: dizziness history of each observation (factorial: yes/no)
- Medication: use of medication of each observation (factorial: yes/no)
- Therapy.Sessions..per.month.: number of therapy sessions of each observation per month (numerical: 1 session)
- Recent.Major.Life.Event: self-reported major life event happened to each observation (factorial: yes/no)
- Diet.Quality..1.10.: diet quality of each observation (factorial: 1-10)
- Anxiety.Level..1.10.: self-reported level of social anxiety of each observation (factorial: 1-10)
Source of the data
https://www.kaggle.com/datasets/natezhang123/social-anxiety-dataset (Kaggle, 2025)
Data manipulation
Despite containing factorial values, the raw dataset does not include any. The main objective of data manipulation is to convert those numeric values into factorial for easier analysis. Furthermore, the column names are renamed for convenience purposes, removing dots from their names. Note that the descriptive statistics only shows the caffeine intake grouped by anxiety level for the purpose of this analysis.
The descriptive statistics shine lights on the structure of the data. For example, the max and min shows the maximal and the minimal value of the group. On the other hand, range shows the difference between max and min, giving an idea as to how spread out the data may be. Median is the number at which \(\leq 50\%\) of the observation is equal or smaller than the value whereas the rest is larger than that value.
From the output, it can be seen that, whilst the max stays at the maximum value, namely around 599 mg caffeine/day, the minimum increases as the anxiety level goes up. This might suggest difference in the distribution of caffeine intake for people reporting different anxiety level. Intuitively, this also reflects in range. As the anxiety level goes up, the range gets smaller, suggesting a narrowing in their distribution. On the other hand, as the range gets smaller, the distribution location changes. This results in a larger median due to distribution location adjusting to the right in the histogram or upwards in the boxplot.
library(psych)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# renaming columns
df <- df %>%
rename(SleepHours = Sleep.Hours) %>%
rename(AnxietyLevel = Anxiety.Level..1.10.) %>%
rename(CaffeineIntake = Caffeine.Intake..mg.day.) %>%
rename(FamilyHistoryOfAnxiety = Family.History.of.Anxiety) %>%
rename(MajorLifeEvent = Recent.Major.Life.Event) %>%
rename(StressLevel = Stress.Level..1.10.) %>%
rename(SweatingLevel = Sweating.Level..1.5.) %>%
rename(DietQuality = Diet.Quality..1.10.) %>%
rename(PhysicalActivityHoursPerWeek = Physical.Activity..hrs.week.) %>%
rename(AlcoholConsumption = Alcohol.Consumption..drinks.week.) %>%
rename(HeartRate = Heart.Rate..bpm.) %>%
rename(BreathingRate = Breathing.Rate..breaths.min.) %>%
rename(TherapySession = Therapy.Sessions..per.month.)
# factorisation of variables
df$Smoking <- factor(df$Smoking, levels = c("Yes", "No"), labels = c("Yes", "No"))
df$FamilyHistoryOfAnxiety <- factor(df$FamilyHistoryOfAnxiety, levels = c("Yes", "No"), labels = c("Yes", "No"))
df$Dizziness <- factor(df$Dizziness, levels = c("Yes", "No"), labels = c("Yes", "No"))
df$Medication <- factor(df$Medication, levels = c("Yes", "No"), labels = c("Yes", "No"))
df$MajorLifeEvent <- factor(df$MajorLifeEvent, levels = c("Yes", "No"), labels = c("Yes", "No"))
df$Occupation <- factor(df$Occupation,
levels = c("Artist", "Nurse", "Other", "Scientist", "Lawyer", "Teacher", "Doctor", "Musician", "Student", "Engineer", "Freelancer", "Chef", "Athlete"),
labels = c("Artist", "Nurse", "Other", "Scientist", "Lawyer", "Teacher", "Doctor", "Musician", "Student", "Engineer", "Freelancer", "Chef", "Athlete"))
df$AnxietyLevel <- as.factor(df$AnxietyLevel)
df$StressLevel <- as.factor(df$StressLevel)
df$SweatingLevel <- as.factor(df$SweatingLevel)
df$DietQuality <- as.factor(df$DietQuality)
# descriptive statistics
#describeBy(df[, -c(1, 2, 3)], group = df$AnxietyLevel) # descriptive statistics on all but demographic variables
describeBy(df[, c(6)], group = df$AnxietyLevel) # descriptive statistics on caffeine intake##
## Descriptive statistics by group
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 1039 239.86 125.87 228 231.01 127.5 2 599 597 0.59 -0.05 3.9
## ------------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 1756 258.57 135.2 250.5 249.28 141.59 1 598 597 0.52 -0.3
## se
## X1 3.23
## ------------------------------------------------------------
## group: 3
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 2407 257.86 136.56 244 248.55 140.85 0 599 599 0.53 -0.33
## se
## X1 2.78
## ------------------------------------------------------------
## group: 4
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 2416 275.47 141.05 263 268.8 151.23 0 599 599 0.37 -0.59
## se
## X1 2.87
## ------------------------------------------------------------
## group: 5
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 1629 284.49 139.3 272 278.14 143.81 2 599 597 0.36 -0.52
## se
## X1 3.45
## ------------------------------------------------------------
## group: 6
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 616 322.75 146.73 294 319.58 147.52 9 599 590 0.25 -0.88
## se
## X1 5.91
## ------------------------------------------------------------
## group: 7
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 123 332.76 149.59 313 334.27 154.19 16 599 583 0.01 -0.82
## se
## X1 13.49
## ------------------------------------------------------------
## group: 8
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 363 444.39 92.05 448 444.72 115.64 179 598 419 -0.1 -0.99 4.83
## ------------------------------------------------------------
## group: 9
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 329 449.4 88.14 440 449.06 108.23 294 598 304 0.08 -1.23 4.86
## ------------------------------------------------------------
## group: 10
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 322 450.94 89.61 459 452.26 117.13 300 599 299 -0.13 -1.31
## se
## X1 4.99Data visualisation
Firstly, the histogram shows the distribution of caffeine consumption by occupation. This allows a quick overview of the distribution. Whilst one can visually infer the normality of each distribution, the empirical result will be examined through Shapiro-Wilk Test.
#install.packages("ggpubr")
library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.2.3## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.2.3##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary(paletteer)## Warning: package 'paletteer' was built under R version 4.2.3gghistogram(df,
x = "CaffeineIntake",
color = "AnxietyLevel",
fill = "AnxietyLevel",
palette = paletteer_d("ggthemes::Classic_Cyclic"),
bins = 100,
title = "Histogram of Daily Caffeine Intake by Occupation",
xlab = "Caffeine Intake (mg/day)",
ylab = "Frequency",
facet.by = "AnxietyLevel") +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(size = 10, hjust = 1),
axis.text.y = element_text(size = 10, hjust = 1),
strip.background = element_rect(fill = "transparent", color = NA),
panel.border = element_blank(),
legend.position = "none",
axis.line.y = element_line(color = "black"))
Secondly, the boxplot below gives an faint overview of whether the homogeneity of variance exists in each group. From the outliers and their spreads, one can detect whether there is a violation of the homoscedasticity. More specifically, if the spread of each group differs from one another, the variance may not be homogeneous. Empirically, tests such as Levene’s test is suitable to statistically confirm this assumption.
ggboxplot(df,
x = "AnxietyLevel",
y = "CaffeineIntake",
color = "AnxietyLevel",
palette = paletteer_d("ggthemes::Classic_Cyclic"),
title = "Boxplot of Daily Caffeine Intake by Occupation",
xlab = "Anxiety Level",
ylab = "Caffeine Intake (mg/day)") +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(size = 12, hjust = 1),
axis.text.y = element_text(size = 10, hjust = 1),
legend.position = "none")
Assumption testing
Shapiro-Wilk test aims to test the normality of the data. Here, the data are grouped by “AnxietyLevel”. Teh two hypotheses are as followed:
\(H_0\): data is normally distributed
\(H_1\): data is not normally distributed
Based on given the p-values at each given anxiety level, the test suggests the data deviates from normality at 95% confidence level. Going back to the histograms above, the fact that none of the graphs follows the shape of a normal distribution can be statistically confirmed by the results below.
library(dplyr)
library(rstatix)## Warning: package 'rstatix' was built under R version 4.2.3##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filtershapiro_results <- df %>%
group_by(AnxietyLevel) %>%
shapiro_test(CaffeineIntake) %>%
mutate(Signif.Code = case_when(
p < 0.001 ~ "***",
p < 0.01 ~ "**",
p < 0.05 ~ "*",
p < 0.1 ~ ".",
TRUE ~ "ns" # not significant
))
print(shapiro_results)## # A tibble: 10 × 5
## AnxietyLevel variable statistic p Signif.Code
## <fct> <chr> <dbl> <dbl> <chr>
## 1 1 CaffeineIntake 0.969 4.22e-14 ***
## 2 2 CaffeineIntake 0.969 5.72e-19 ***
## 3 3 CaffeineIntake 0.968 1.12e-22 ***
## 4 4 CaffeineIntake 0.974 1.84e-20 ***
## 5 5 CaffeineIntake 0.976 8.10e-16 ***
## 6 6 CaffeineIntake 0.963 2.07e-11 ***
## 7 7 CaffeineIntake 0.973 1.59e- 2 *
## 8 8 CaffeineIntake 0.960 2.37e- 8 ***
## 9 9 CaffeineIntake 0.949 3.23e- 9 ***
## 10 10 CaffeineIntake 0.939 3.36e-10 ***Levene test aims to examine the homoscedasticity of the variable. The test assumes:
\(H_0\): the variance across all occupation is the same
\(H_1\): at least one variance is different from others
Having examined the output, Levene test concludes the p-value as \(<2.2*10^{-16}\). The p-value is significantly smaller than that required at 99% significant level. Based on the output, one can infer that the violation of homoscedasticity is statistically significant. As such, it might be more suitable to employ Welch t-test instead of ANOVA given the former is less affected by heteroscedasticity.
library(car)## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:psych':
##
## logitleveneTest(df$CaffeineIntake, group = df$AnxietyLevel)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 9 20.604 < 2.2e-16 ***
## 10990
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Parametric test (Welch test)
Welch test is a parametric test that compares the means of multiple groups. Unlike standard ANOVA, Welch test can perform more robustly when heteroscedasticity is found.
The test assumes the following hypotheses:
\(H_0\): all group means are equal
\(H_1\): at least one group mean is different
Based on the result, the null hypothesis can be rejected and it can be determined that, at 95% confidence level, at least one group mean is different.
#install.packages("onewaytests")
library(onewaytests)## Warning: package 'onewaytests' was built under R version 4.2.3##
## Attaching package: 'onewaytests'## The following object is masked from 'package:psych':
##
## describewelch_results <- welch.test(CaffeineIntake ~ AnxietyLevel,
data = df)##
## Welch's Heteroscedastic F Test (alpha = 0.05)
## -------------------------------------------------------------
## data : CaffeineIntake and AnxietyLevel
##
## statistic : 386.8586
## num df : 9
## denom df : 1685.508
## p.value : 0
##
## Result : Difference is statistically significant.
## -------------------------------------------------------------Pairwise T-test pairs each possible group and tests whether there is a difference between their means. Since Welch test only tests whether there is difference in means at all within the data, pairwise test allows us examine closer. It can be determined that the caffeine intake between people who have responded with 2 and 3 anxiety level does not have significant difference in their means. The fact that they have chosen to report different anxiety level may be motivated by other factors.
Below is a list of pairs where null-hypothesis cannot be rejected at \(P < 0.05\) after Bonferroni adjustment.
- (2, 3)
- (4, 5)
- (6, 7)
- (8, 9)
- (8, 10)
- (9, 10)
pairwise.t.test(x = df$CaffeineIntake,
g = df$AnxietyLevel,
p.adjust.method = "bonf")##
## Pairwise comparisons using t tests with pooled SD
##
## data: df$CaffeineIntake and df$AnxietyLevel
##
## 1 2 3 4 5 6 7 8 9
## 2 0.01629 - - - - - - - -
## 3 0.01337 1.00000 - - - - - - -
## 4 3.7e-11 0.00259 0.00023 - - - - - -
## 5 2.4e-15 8.5e-07 2.7e-08 1.00000 - - - - -
## 6 < 2e-16 < 2e-16 < 2e-16 2.6e-13 7.2e-08 - - - -
## 7 1.7e-11 1.3e-07 6.8e-08 0.00017 0.00525 1.00000 - - -
## 8 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 6.8e-14 - -
## 9 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 8.7e-15 1.00000 -
## 10 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 4.3e-15 1.00000 1.00000
##
## P value adjustment method: bonferroniBelow is the graphical presentation of Welch Test.
library(dplyr)
library(rstatix)
pwc_pw <- df %>%
pairwise_t_test(CaffeineIntake ~ AnxietyLevel,
paired = FALSE,
p.adjust.method = "bonferroni")
welch_results <- welch_anova_test(CaffeineIntake ~ AnxietyLevel,
data = df)
pwc_pw <- pwc_pw %>%
add_y_position(fun = "median", step.increase = 0.35)
ggboxplot(df, x = "AnxietyLevel", y = "CaffeineIntake",
color = "AnxietyLevel", palette = paletteer_d("ggthemes::Classic_Cyclic")) +
coord_cartesian(ylim = c(0, 600)) + # Focus on relevant range
stat_summary(fun = mean, geom = "point", shape = 16, size = 3,
aes(group = AnxietyLevel), color = "darkred") +
stat_summary(fun = mean, geom = "text", colour = "darkred", size = 4,
vjust = -1, hjust = 0, aes(label = round(after_stat(y), 2), group = AnxietyLevel)) +
theme(axis.text.x = element_text(hjust = 1, size = 12),
legend.position = "none") + # Improve x-axis readability
labs(subtitle = get_test_label(welch_results, detailed = TRUE),
caption = get_pwc_label(pwc_pw))
Non-parametric test (Kruskal-Wallis Test)
As shown above, the normality of the data appears to be violated. It might be better then, therefore, to conduct the non-parametric test, namely the Kruskal-Wallis Rank Sum Test. The test aims to compare the distribution location of the variables. The hypotheses are as followed:
\(H_0\): the distribution location of caffeine intake is the same across all groups
\(H_1\): at least one of the distribution location of caffeine intake is different
Based on the output, the null-hypothesis is rejected at \(p < 0.001\).
kruskal.test(CaffeineIntake ~ AnxietyLevel,
data = df)##
## Kruskal-Wallis rank sum test
##
## data: CaffeineIntake by AnxietyLevel
## Kruskal-Wallis chi-squared = 1550.9, df = 9, p-value < 2.2e-16Subsequently, the output below shows the effect size, namely how different are the distribution locations. The output suggests that the difference is large.
kruskal_effsize(CaffeineIntake ~ AnxietyLevel,
data = df)## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 CaffeineIntake 11000 0.140 eta2[H] largeAs steps above, the testing can be followed by a pair-wise approach. Wilcoxon rank-sum test is chosen, since it’s suitable for multiple comparisons of non-parametric test.
Below is a list of paired groups of which the null-hypothesis cannot be rejected at \(p < 0.05\).
- (1, 3)
- (2, 3)
- (4, 5)
- (6, 7)
- (8, 9)
- (8, 10)
- (9, 10)
Different from the output by pairwise t-test, given the normality assumption, the output from Wilcoxon rank sum test shows less distribution location to be deviating from one another. One of the main differences is the fact that Wilcoxon test takes median as its assumption, neglecting the distribution of the data.
group_nonpar <- wilcox_test(CaffeineIntake ~ AnxietyLevel,
paired = FALSE,
p.adjust.method = "bonferroni",
data = df)
group_nonpar ## # A tibble: 45 × 9
## .y. group1 group2 n1 n2 statis…¹ p p.adj p.adj…²
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 CaffeineIntake 1 2 1039 1756 843578 8.68e- 4 3.9 e- 2 *
## 2 CaffeineIntake 1 3 1039 2407 1163954. 1 e- 3 5.6 e- 2 ns
## 3 CaffeineIntake 1 4 1039 2416 1074619 1.91e- 11 8.6 e- 10 ****
## 4 CaffeineIntake 1 5 1039 1629 687856 3.23e- 16 1.45e- 14 ****
## 5 CaffeineIntake 1 6 1039 616 216708 4.17e- 28 1.88e- 26 ****
## 6 CaffeineIntake 1 7 1039 123 40412. 2.5 e- 11 1.12e- 9 ****
## 7 CaffeineIntake 1 8 1039 363 38847 1.39e-112 6.26e-111 ****
## 8 CaffeineIntake 1 9 1039 329 32080. 1.67e-109 7.52e-108 ****
## 9 CaffeineIntake 1 10 1039 322 32030. 8.98e-107 4.04e-105 ****
## 10 CaffeineIntake 2 3 1756 2407 2124961 7.62e- 1 1 e+ 0 ns
## # … with 35 more rows, and abbreviated variable names ¹statistic, ²p.adj.signifBelow is the visualisation of Kruskal-Wallis test.
pwc_wt <- df %>%
wilcox_test(CaffeineIntake ~ AnxietyLevel,
paired = FALSE,
p.adjust.method = "bonferroni")
Kruskal_results <- kruskal_test(CaffeineIntake ~ AnxietyLevel,
data = df)
library(rstatix)
pwc_wt <- pwc_wt %>%
add_y_position(fun = "max", step.increase = 0.20)
library(ggpubr)
ggboxplot(df, x = "AnxietyLevel", y = "CaffeineIntake", color = "AnxietyLevel", palette = paletteer_d("ggthemes::Classic_Cyclic")) +
coord_cartesian(ylim = c(0, 600)) +
stat_summary(fun = median, geom = "point", shape = 16, size = 3,
aes(group = AnxietyLevel), color = "darkred") +
stat_summary(fun = median, geom = "text", colour = "darkred", size = 4,
vjust = -1, hjust = 0, aes(label = round(after_stat(y), 2), group = AnxietyLevel)) +
theme(axis.text.x = element_text(hjust = 1, size = 12),
legend.position = "none") + # Improve x-axis readability
labs(subtitle = get_test_label(Kruskal_results, detailed = TRUE),
caption = get_pwc_label(pwc_wt))
Conclusion
Based on the results, it can be concluded that Kruskal-Wallis test might have captured the nature of the data slightly better than Welch test. From the assumption testing, the data is concluded to be in violation of both normality and homoscedasticity, compromising the statistical power of a parametric test. The Kruskal-Wallis detects less variation in caffeine intake across different pairs of anxiety level, which may suggest that caffeine intake is not the only variable that influences respondent’s anxiety level. However, it may be worth noting that the size of each group is quite different from one another. The paired test may have captured the slight changes in medians of larger samples and concluded as such. As such, it might be worth to standardise the sample size of each group and repeat the analysis. In conclusion, the non-parametric test may be more suitable in capturing the variation in distribution location, suggesting fewer differences in caffeine intake for certain pairs of anxiety levels.