Internet use, religiosity and happiness in Netherlands

Introduction

In this project I have decided to explore a variable which can act as a mediator between being religious (belonging to a certain religion, visiting religious events, etc.) and being happy. This variable is netustm (Internet use, how much time on typical day, in minutes). This decision was inspired by this research. The results of it were as follows: “those who reported spending less time on the internet, less time expressing emotions, and more time checking facts, scored higher on measures of happiness”. I also have decided to check my own hypothesis on whether discriminated people spent more time on the Internet. Behind this idea lies the hypothesis that people who belong to a discriminated groups are more likely to seek community online because of safety concerns and feeling of isolation.

Data preparation and other technical stuff

Activate packages

library(kableExtra)
library(readr)
library(dplyr)
library(ggplot2)
library(tidyr)
library(moments)
library(perm)
library(effectsize)
library(car)
library(psych)
library(stats)
library(sjstats)
library(pwr)
library(DescTools)
library(rcompanion)

Import the data set

data <- read.csv("C:/Users/Asus/Downloads/Data.csv")

Select target variables

data <- select(data, c(cntry, rlgblg, netustm, dscrgrp, netusoft, rlgatnd))

Create a data set with observations from Netherlands only

ESS_NL <- data %>% filter(cntry == "NL")

Eliminate unwanted values while preserving the initial amount of observations, assign labels, set a correct type of data for each variable

ESS_NL$rlgblg <- ifelse (ESS_NL$rlgblg > 2, NA, ifelse(ESS_NL$rlgblg > 1, "No", "Yes"))
ESS_NL$netusoft <- ifelse (ESS_NL$netusoft <= 2, "Rarely", 
                         ifelse (ESS_NL$netusoft < 4, "Moderately", 
                         ifelse (ESS_NL$netusoft > 5, NA, "Often")))
ESS_NL$rlgatnd <- ifelse (ESS_NL$rlgatnd <= 3, "Often", 
                         ifelse (ESS_NL$rlgatnd < 5, "Moderately", 
                          ifelse (ESS_NL$rlgatnd > 7, NA, "Rarely")))
ESS_NL$netustm <- ifelse(ESS_NL$netustm > 1440, NA, ESS_NL$netustm)
ESS_NL$dscrgrp <- ifelse (ESS_NL$dscrgrp > 2, NA, ifelse(ESS_NL$dscrgrp > 1, "No", "Yes"))
ESS_NL$rlgblg <- factor(ESS_NL$rlgblg, ordered = FALSE)
ESS_NL$netusoft <- factor(ESS_NL$netusoft, ordered = FALSE)
ESS_NL$dscrgrp <- factor(ESS_NL$dscrgrp, ordered = FALSE)
ESS_NL$rlgatnd <- factor(ESS_NL$rlgatnd, ordered = FALSE)

Chi-square

Variables rlgblg and eisced are both categorical, which makes them suitable for chi-square test. Our alternative hypothesis is that belonging to a certain religion and level of education are connected, therefore our null hypothesis is that there is no connection between those variables.

Create a data set

ESS_chi <- ESS_NL %>% 
  select(rlgblg, netusoft)

Create a cross-table

cross <- table(ESS_chi$netusoft, ESS_chi$rlgblg)
kable(cross)%>% 
  kable_styling(bootstrap_options=c("bordered", "responsive","striped"), full_width = FALSE)

	No	Yes
Moderately	11	14
Often	1030	357
Rarely	40	13

Create a stacked bar plot with variables used

ggplot(data = ESS_chi[!(is.na(ESS_chi$rlgblg)),], aes(x = netusoft, fill=rlgblg)) +
  geom_bar()+
  xlim("Rarely", "Moderately", "Often") +
  coord_flip()+
  xlab("Internet use, how often") + 
  ylab("Belonging to a religion") +
  ggtitle("Belonging to a religion and the Internet use") +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_fill_brewer(palette = "PuBuGn", breaks=c('Yes', 'No'))+
  theme_minimal()

Before performing a chi-square test we should check if the assumptions for it are met:

“The data in the cells should be frequencies, or counts of cases rather than percentages or some other transformation of the data.” - check
“The levels (or categories) of the variables are mutually exclusive.” - check
“Each subject may contribute data to one and only one cell in the χ2.” - check
“The study groups must be independent.” - check

Since all of our assumptions are met we can now proceed with the test:

Perform a chi-square test

chisq <- chisq.test(ESS_chi$rlgblg, ESS_chi$netusoft)
chisq

## 
##  Pearson's Chi-squared test
## 
## data:  ESS_chi$rlgblg and ESS_chi$netusoft
## X-squared = 11.708, df = 2, p-value = 0.002869

P-value is <0,05, therefore we can reject the null hypothesis.

Now we should analyse standardized residuals:

Analyse standardized residuals

kable(chisq$stdres)%>% 
  kable_styling(bootstrap_options=c("bordered", "responsive","striped"), full_width = FALSE)

	Moderately	Often	Rarely
No	-3.415968	1.73445	0.2838315
Yes	3.415968	-1.73445	-0.2838315

There are categories (religious who use Internet moderately and non-religious who use the Internet moderately) which have residuals >2 or <-2, which allows us to say that those categories are influencing chi-square the most.

Conclusion

We can say that there is a connection between person’s belonging to a certain religion and the frequency of the Internet use.

T-test

Variable dscrgrp is categorical (it is going to be our independent variable) and variable netustm is numerical (dependent variable), which makes them suitable for t-test. Our alternative hypothesis is members of discriminated groups on average spent more time on the Internet than those who do not belong to any discriminated groups, therefore our null hypothesis is that there is no significant difference between those two groups.

Create a data set

ESS_tt <- ESS_NL %>% 
  select(dscrgrp, netustm)

Create a box plot with a visualization of distribution in groups

ggplot() +
  geom_boxplot(data = ESS_tt[!(is.na(ESS_tt$dscrgrp)),], aes(x = dscrgrp, y = netustm), fill="#a6bddb", col="#3690c0", alpha = 0.5) +
  xlab("Belongs to a discriminated group") + 
  ylab("Amount of time spent on the Internet (minutes)") +
  ggtitle("Belonging to a discriminated group and time spent on the Internet") +
  theme_minimal()

We can already see that members of discriminated groups have higher mean.

Before performing a t-test we should check if the assumptions for it are met:

“Observations are independent” - check
The distribution is normal/close to normal - in order to check this we have to explore certain markers of normality:

Skewness and kurtosis

skewness(ESS_tt$netustm, na.rm = TRUE)

## [1] 0.9691282

kurtosis(ESS_tt$netustm, na.rm = TRUE)

## [1] 4.227383

Skewness is positive, it indicates that the distribution is right-skewed; Kurtosis is bigger than 2, it indicates that the distribution is very peaked and has more values in the tail compared to normal distribution.

Histogram

ggplot(ESS_tt[!(is.na(ESS_tt$dscrgrp)),], aes(x = netustm, fill = dscrgrp)) +
  geom_histogram(aes(y=..density..), position = "identity", alpha = 0.7, binwidth = 50) +
  geom_density(col = "black", fill = "white", alpha = 0.1) +
  xlab("Amount of time spent on the Internet (minutes)") + 
  ylab("Density") +
  ggtitle("Belonging to a discriminated group and time spent on the Internet") +
  scale_fill_brewer(palette = "PuBuGn", breaks=c('Yes', 'No'))+
  theme_minimal()

The distribution is quite abnormal, both for discriminated and non-discriminated. However, amount of observations is >300, so we can neglect the abnormality of distribution because of the central limit theorem (we can also use permutation test and Wilcoxon test for better results).

“Equal variances across two groups” - got to perform additional tests as well:

Levene’s test

leveneTest(ESS_tt$netustm ~ ESS_tt$dscrgrp)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    1  2.4388 0.1186
##       1375

P-value in this case is > 0.05, so we cannot reject H0. It means that variances are not distributed equally, therefore we have to perform Welch’s t-test instead of a regular one.

Perform a Welch’s t-test

t.test(ESS_tt$netustm ~ ESS_tt$dscrgrp, var.equal = F)

## 
##  Welch Two Sample t-test
## 
## data:  ESS_tt$netustm by ESS_tt$dscrgrp
## t = -4.2143, df = 220.89, p-value = 3.652e-05
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
##  -113.96735  -41.34017
## sample estimates:
##  mean in group No mean in group Yes 
##          292.1474          369.8011

P-value is < 0.05, therefore we can reject the null hypothesis.

Since the results are statistically significant we should report the effect size:

Calculate the effect size

cohens_d(ESS_tt$netustm ~ ESS_tt$dscrgrp)

## Cohen's d |         95% CI
## --------------------------
## -0.36     | [-0.52, -0.20]
## 
## - Estimated using pooled SD.

Effect size of -0.36 can be labeled as small.

Now we should double-check our results with a non-parametric test (I have chosen a permutation test):

Perform a non-parametric test

permTS(netustm ~ dscrgrp, data = ESS_tt)

## 
##  Permutation Test using Asymptotic Approximation
## 
## data:  netustm by dscrgrp
## Z = -4.4439, p-value = 8.836e-06
## alternative hypothesis: true mean dscrgrp=No - mean dscrgrp=Yes is not equal to 0
## sample estimates:
## mean dscrgrp=No - mean dscrgrp=Yes 
##                          -77.65376

P-value from permutations test is < 0.05 as well, which also allows as to reject H0.

Conclusion

Delta between means (mean “No”- mean “Yes”) is negative and results of t-test and permutations test are statistically significant, therefore we can say that members of discriminated groups on average spent more time on the Internet than those who do not belong to any discriminated groups (with caution, however, since the effect size is marked as small).

ANOVA

Variable rlgatnd is categorical with three levels (independent variable) and variable netustm is numerical (dependent variable), which makes them suitable for ANOVA test. Our alternative hypothesis is that there is a significant difference in the average time spent on the Internet between at least two groups attending religious events, therefore our null hypothesis is that there is no difference at all.

Create a data set with religious people only

ESS_anova <- ESS_NL %>% filter(rlgblg == "Yes")
ESS_anova <- select(ESS_anova, c(rlgatnd, netustm))

Create a table with a statistical description of groups

description <- ESS_anova[!is.na(ESS_anova$rlgatnd),] %>%
  group_by(rlgatnd) %>%
  summarise(N = n(), Mean = mean(netustm, na.rm = T), SD = sd(netustm, na.rm = TRUE))
kable(description)%>% 
  kable_styling(bootstrap_options=c("bordered", "responsive","striped"), full_width = FALSE)

rlgatnd	N	Mean	SD
Moderately	56	187.7885	154.3056
Often	116	233.6226	198.9144
Rarely	211	293.3163	240.6364

As well as in the previous case, we can see that there is a difference between group means.

Create a box plot with a visualization of distribution in groups

ggplot() +
  geom_boxplot(data = ESS_anova[!is.na(ESS_anova$rlgatnd),], aes(x = rlgatnd, y = netustm), fill="#a6bddb", col="#3690c0") +
  xlab("Attendance of religious events") + 
  ylab("Time spent on the Internet") +
  ggtitle("Use of the Internet and attendance of religious events")+
  theme(plot.title = element_text(hjust = 0.5))

Before performing an ANOVA test we should check if the assumptions for it are met:

“Independence of observations” - check
“Normal distribution of residuals” - going to check it after the performance of a f-test
“Equal variances” - in order to check this we should once again run a Levene’s test:

leveneTest(ESS_anova$netustm ~ ESS_anova$rlgatnd)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value   Pr(>F)   
## group   2  5.4497 0.004669 **
##       351                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

P-value is < 0.05, therefore we cannot reject H0 for this test. It means that there is a significant difference between variables in groups. We should adjust our f-test accordingly.

Perform a f-test

oneway.test(ESS_anova$netustm ~ ESS_anova$rlgatnd, var.equal = F)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  ESS_anova$netustm and ESS_anova$rlgatnd
## F = 7.6624, num df = 2.00, denom df = 159.46, p-value = 0.0006647

P-value is < 0.05 therefore we can reject H0. It means that there is a significant difference between means of at least two groups. Now we should check if residuals are normally distributed:

Check the normality of residuals’ distribution

one.way.anova <- aov(ESS_anova$netustm ~ ESS_anova$rlgatnd) #anova was not working
summary(one.way.anova)

##                    Df   Sum Sq Mean Sq F value  Pr(>F)   
## ESS_anova$rlgatnd   2   562460  281230   5.925 0.00295 **
## Residuals         351 16660492   47466                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 30 пропущенных наблюдений удалены

plot(one.way.anova, 2)

anova.res <- residuals(object = one.way.anova)
describe(anova.res) # skew < 2, kurtosis > 2

##    vars   n mean     sd median trimmed    mad     min     max range skew
## X1    1 354    0 217.25 -53.62   -25.7 177.91 -278.32 1146.68  1425  1.5
##    kurtosis    se
## X1     3.75 11.55

shapiro.test(x = anova.res)

## 
##  Shapiro-Wilk normality test
## 
## data:  anova.res
## W = 0.87446, p-value = 2.266e-16

On the plot we can see that distribution is somewhat close to normal, however Shapiro-Wilk test says otherwise (p<0,05), and kurtosis is >2, which also indicates abnormality.

Perform a post-hoc test (Bonferroni correction is used since the variances are not equal)

pairwise.t.test(ESS_anova$netustm, ESS_anova$rlgatnd, 
                p.adjust.method = "bonferroni", pool.sd = T)

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  ESS_anova$netustm and ESS_anova$rlgatnd 
## 
##        Moderately Often 
## Often  0.6446     -     
## Rarely 0.0062     0.0710
## 
## P value adjustment method: bonferroni

There is one p-value < 0.05, it is for the difference between rarely and moderately, so we can say that the only significant difference is between people who attend religious events rarely and those who attend moderately.

Calculate effect size

stats1 <- anova_stats(one.way.anova)

Effect size of 0.027 can be labeled as small.

Perform a non-parametric test

kruskal.test(netustm ~ rlgatnd, data = ESS_anova)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  netustm by rlgatnd
## Kruskal-Wallis chi-squared = 10.635, df = 2, p-value = 0.004904

P-value is < 0.05 therefore we can reject H0. It confirms results provided by an ANOVA test.

Perform a post-hoc test

DunnTest(netustm ~ rlgatnd, data = ESS_anova)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                   mean.rank.diff   pval    
## Often-Moderately        20.43015 0.2366    
## Rarely-Moderately       46.83379 0.0097 ** 
## Rarely-Often            26.40364 0.0634 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The results are replicating those which were provided by a paired t-test.

Calculate effect size

epsilonSquared(x = ESS_anova$netustm, g = ESS_anova$rlgatnd)

## epsilon.squared 
##          0.0278

Effect size is very close to one reported in step 7 and can be marked as small.

Conclusion

The results of an ANOVA test an a Kruskall-Wallis test are statistically significant, therefore we can say that there is a significant difference in the average time spent on the Internet between at least two groups attending religious events (with caution, however, since the effect size is marked as small).

Final conclusion

All of the hypotheses were confirmed, however most of them reportedly have a small effect size, which lefts me a little suspicious about the results. I will develop used ideas further in our next project.

Internet use, religiosity and happiness in Netherlands

Daria Shevchenko

2024-03-04

Introduction

Data preparation and other technical stuff

Chi-square

Conclusion

T-test

Conclusion

ANOVA

Conclusion

Final conclusion