Main analysis
1. Time trends of life satisfaction
Extract means & proportions from each calendar year (intercept-only models)
# TAKES A FEW MINS TO RUN! Output saved below
dvlist <- c("h.school.c", "h.schoolwork.c", "h.friends.c","h.fam.c", "h.appear.c", "h.life.c" )
dvlist2 <- c("unhappy.school", "unhappy.schoolwk", "unhappy.friends","unhappy.fam", "unhappy.appear", "unhappy.life" )
## FULL SAMPLE #####
means.fullsample <-c()
prop.fullsample <- c()
# FULL
for (year.x in sort(unique(usoc$int.year.sched)) ){ # loops through interview years (in order!)
xs <- usoc[usoc$int.year.sched==year.x , ] # isolate data for one year
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled, # note if wgt=0, the person is dropped from estimates. But they may be retained for cluster and stratum calculations.
data=xs,
nest=TRUE)
options(survey.lonely.psu = "adjust") # tells the survey package to handle singleton PSUs by scaling their variance estimates
# run svyglm with intercept-only to estimate population-weighted *mean* for each domain this wave
dom.avg.m <- lapply(dvlist, function(x) {
svyglm(eval(substitute(i ~ 1 ,
list(i = as.name(x)))),
design = usoc.svy)})
# Warning :observations with zero weight not used for calculating dispersion
# separately calculate the weighted population variance-- how much does satisfaction vary (across individuals) in the population
dom.var <- lapply(dvlist, function(x) {
svyvar(eval(substitute(~ i ,
list(i = as.name(x)))),
design = usoc.svy)})
# extract mean and se from the intercept-only model
x <- c()
for (domain in 1:6){
x <- as.data.frame(rbind(x,
data.frame(year.x=year.x,
domain=dvlist[domain],
mean = summary(dom.avg.m[[domain]])[["coefficients"]][1], # mean
SE = summary(dom.avg.m[[domain]])[["coefficients"]][2], # this is the sampling error of the mean.
# a measure of precision of the mean.
popVar = as.data.frame(dom.var[[domain]])[,1], # *this* is the population variance, a measure of *spread* and inequality
popVar.SE = as.data.frame(dom.var[[domain]])[,2],
N = length(dom.avg.m[[1]][["effects"]])
)))
}
# find weighted proportion of unhappy faces in this year (with 95% CI)
prop <- lapply(dvlist2, function(x) {
svyciprop(eval(substitute(i ~ 1 ,
list(i = as.name(x)))),
design = usoc.svy)})
y <- c()
for (domain in 1:6){
y <- as.data.frame(rbind(y,
data.frame( cbind(year.x,
names(prop[[domain]]),
prop[[domain]],
confint(prop[[domain]]) ) )
))
}
rownames(y) <- c()
names(y) <- c("year.x", "domain", "prop", "ci.low", "ci.hi")
y$year.x <- as.numeric(y$year.x)
y$prop <- as.numeric(y$prop)
y$ci.low <- as.numeric(y$ci.low)
y$ci.hi <- as.numeric(y$ci.hi)
means.fullsample <- rbind(means.fullsample, x)
prop.fullsample <- rbind(prop.fullsample, y)
} # next wave
# centre time in years such that 0=2009 (for linear trend across full period )
means.fullsample$year.n <- as.numeric(as.factor(means.fullsample$year.x))-1
# break years into 3 windows, each with 5 years in them: 2010-2014; 2015-2019; 2020-2024 (i.e. remove 2009)
means.fullsample$windows <- cut(means.fullsample$year.x, breaks=c(2009, 2014, 2019, 2025), labels=c(1, 2, 3))
# repeat for the proportions dataframe
prop.fullsample$year.n <- as.numeric(as.factor(prop.fullsample$year.x))-1
prop.fullsample$windows <- cut(prop.fullsample$year.x, breaks=c(2009, 2014, 2019, 2025), labels=c(1, 2, 3))
# extract weight means for each subgroup (males, females, age 10-12, 13-15, low, middle & high income)
# BY SEX #####
means.bysex <- c()
props.bysex <- c()
for (sex in c("Female", "Male")) {
all.res <- c() # holds all results for one group (e.g. females)
for (year.x in sort(unique(usoc$int.year.sched)) ){
xs <- usoc[usoc$int.year.sched==year.x,] # isolate data for this YEAR
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
tmp <- usoc.svy$variables$sex_dv == sex # limit survey data to group
usoc.svg.group <- subset(usoc.svy, tmp)
options(survey.lonely.psu = "adjust") # tells the survey package to handle singleton PSUs by scaling their variance estimates
group.glm <- lapply(dvlist, function(x) {
svyglm(eval(substitute(i ~ 1,
list(i = as.name(x)))),
design = usoc.svg.group)})
# extract cross-sectional results for all domains in this YEAR and for this GROUP
xs.res <- c()
for (domain in 1:6){
xs.res <- as.data.frame(rbind(xs.res,
cbind(year.x, sex, dvlist[domain],
t(summary(group.glm[[domain]])[["coefficients"]][1:2]))))
}
# find weighted proportion of unhappy faces in this year (with 95% CI)
prop <- lapply(dvlist2, function(x) {
svyciprop(eval(substitute(i ~ 1 ,
list(i = as.name(x)))),
design = usoc.svg.group)})
y <- c()
for (domain in 1:6){
y <- as.data.frame(rbind(y,
data.frame( cbind(sex,
year.x,
names(prop[[domain]]),
prop[[domain]],
confint(prop[[domain]]) ) )))
}
rownames(y) <- c()
names(y) <- c("year.x","group", "domain", "prop", "ci.low", "ci.hi")
y$year.x <- as.numeric(y$year.x)
y$prop <- as.numeric(y$prop)
y$ci.low <- as.numeric(y$ci.low)
y$ci.hi <- as.numeric(y$ci.hi)
names(xs.res) <- c("year.x", "group", "dvlist", "mean", "se")
all.res <- rbind(all.res, xs.res) # merge results from all years
all.res$mean <- as.numeric(all.res$mean)
all.res$se <- as.numeric(all.res$se)
} # finishes analysis of one year
means.bysex <- rbind(means.bysex, all.res)
props.bysex <- rbind(props.bysex, y)
} # finishes analyses for both groups
means.bysex$year.x <- as.numeric(means.bysex$year.x)
# BY AGE-GROUP #####
means.byage <- c()
props.byage <- c()
for (age in c("10-12", "13-15")) {
all.res <- c() # holds all results for one group (e.g. females)
for (year.x in sort(unique(usoc$int.year.sched)) ){
xs <- usoc[usoc$int.year.sched==year.x,] # isolate data for this YEAR
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
tmp <- usoc.svy$variables$age.group == age # limit survey data to group
usoc.svg.group <- subset(usoc.svy, tmp)
options(survey.lonely.psu = "adjust") # tells the survey package to handle singleton PSUs by scaling their variance estimates
group.glm <- lapply(dvlist, function(x) {
svyglm(eval(substitute(i ~ 1,
list(i = as.name(x)))),
design = usoc.svg.group)})
# extract cross-sectional results for all domains in this YEAR and for this GROUP
xs.res <- c()
for (domain in 1:6){
xs.res <- as.data.frame(rbind(xs.res,
cbind(year.x, age, dvlist[domain],
t(summary(group.glm[[domain]])[["coefficients"]][1:2]))))
}
names(xs.res) <- c("year.x", "group", "dvlist", "mean", "se")
all.res <- rbind(all.res, xs.res) # merge results from all years
all.res$mean <- as.numeric(all.res$mean)
all.res$se <- as.numeric(all.res$se)
# find weighted proportion of unhappy faces in this year (with 95% CI)
prop <- lapply(dvlist2, function(x) {
svyciprop(eval(substitute(i ~ 1 ,
list(i = as.name(x)))),
design = usoc.svg.group)})
y <- c()
for (domain in 1:6){
y <- as.data.frame(rbind(y,
data.frame( cbind(year.x,
age,
names(prop[[domain]]),
prop[[domain]],
confint(prop[[domain]]) ) )))
}
rownames(y) <- c()
names(y) <- c("year.x", "group", "domain", "prop", "ci.low", "ci.hi")
y$year.x <- as.numeric(y$year.x)
y$prop <- as.numeric(y$prop)
y$ci.low <- as.numeric(y$ci.low)
y$ci.hi <- as.numeric(y$ci.hi)
} # finishes analysis for all years
means.byage <- rbind(means.byage, all.res)
props.byage <- rbind(props.byage, y)
} # finishes analyses for both groups
means.byage$year.x <- as.numeric(means.byage$year.x)
# BY SPECIFIC AGE #####
means.byage2 <- c()
props.byage2 <-c()
for (age in c(10,11,12,13,14,15)) {
all.res <- c() # holds all results for one group (e.g. females)
for (year.x in sort(unique(usoc$int.year.sched)) ){
xs <- usoc[usoc$int.year.sched==year.x,] # isolate data for this YEAR
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
tmp <- usoc.svy$variables$dvage.n == age # limit survey data to group
usoc.svg.group <- subset(usoc.svy, tmp)
options(survey.lonely.psu = "adjust") # tells the survey package to handle singleton PSUs by scaling their variance estimates
group.glm <- lapply(dvlist, function(x) {
svyglm(eval(substitute(i ~ 1,
list(i = as.name(x)))),
design = usoc.svg.group)})
# extract cross-sectional results for all domains in this YEAR and for this GROUP
xs.res <- c()
for (domain in 1:6){
xs.res <- as.data.frame(rbind(xs.res,
cbind(year.x, age, dvlist[domain],
t(summary(group.glm[[domain]])[["coefficients"]][1:2]))))
}
names(xs.res) <- c("year.x", "group", "dvlist", "mean", "se")
all.res <- rbind(all.res, xs.res) # merge results from all years
all.res$mean <- as.numeric(all.res$mean)
all.res$se <- as.numeric(all.res$se)
# find weighted proportion of unhappy faces in this year (with 95% CI)
prop <- lapply(dvlist2, function(x) {
svyciprop(eval(substitute(i ~ 1 ,
list(i = as.name(x)))),
design = usoc.svg.group)})
y <- c()
for (domain in 1:6){
y <- as.data.frame(rbind(y,
data.frame( cbind(year.x,
age,
names(prop[[domain]]),
prop[[domain]],
confint(prop[[domain]]) ) )))
}
rownames(y) <- c()
names(y) <- c("year.x","group", "domain", "prop", "ci.low", "ci.hi")
y$year.x <- as.numeric(y$year.x)
y$prop <- as.numeric(y$prop)
y$ci.low <- as.numeric(y$ci.low)
y$ci.hi <- as.numeric(y$ci.hi)
} # finishes analysis for all years
means.byage2 <- rbind(means.byage2, all.res)
props.byage2 <- rbind(props.byage2, y)
} # finishes analyses for both groups
means.byage2$year.x <- as.numeric(means.byage2$year.x)
# BY SPECIFIC AGE (3-group) #####
means.byage3 <- c()
props.byage3 <-c()
for (age in c("10-11", "12-13", "14-15")) {
all.res <- c() # holds all results for one group (e.g. females)
for (year.x in sort(unique(usoc$int.year.sched)) ){
xs <- usoc[usoc$int.year.sched==year.x,] # isolate data for this YEAR
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
tmp <- usoc.svy$variables$age.3groups == age # limit survey data to group
usoc.svg.group <- subset(usoc.svy, tmp)
options(survey.lonely.psu = "adjust") # tells the survey package to handle singleton PSUs by scaling their variance estimates
group.glm <- lapply(dvlist, function(x) {
svyglm(eval(substitute(i ~ 1,
list(i = as.name(x)))),
design = usoc.svg.group)})
# extract cross-sectional results for all domains in this YEAR and for this GROUP
xs.res <- c()
for (domain in 1:6){
xs.res <- as.data.frame(rbind(xs.res,
cbind(year.x, age, dvlist[domain],
t(summary(group.glm[[domain]])[["coefficients"]][1:2]))))
}
names(xs.res) <- c("year.x", "group", "dvlist", "mean", "se")
all.res <- rbind(all.res, xs.res) # merge results from all years
all.res$mean <- as.numeric(all.res$mean)
all.res$se <- as.numeric(all.res$se)
# find weighted proportion of unhappy faces in this year (with 95% CI)
prop <- lapply(dvlist2, function(x) {
svyciprop(eval(substitute(i ~ 1 ,
list(i = as.name(x)))),
design = usoc.svg.group)})
y <- c()
for (domain in 1:6){
y <- as.data.frame(rbind(y,
data.frame( cbind(year.x,
age,
names(prop[[domain]]),
prop[[domain]],
confint(prop[[domain]]) ) )))
}
rownames(y) <- c()
names(y) <- c("year.x","group", "domain", "prop", "ci.low", "ci.hi")
y$year.x <- as.numeric(y$year.x)
y$prop <- as.numeric(y$prop)
y$ci.low <- as.numeric(y$ci.low)
y$ci.hi <- as.numeric(y$ci.hi)
} # finishes analysis for all years
means.byage3 <- rbind(means.byage3, all.res)
props.byage3 <- rbind(props.byage3, y)
} # finishes analyses for both groups
means.byage3$year.x <- as.numeric(means.byage3$year.x)
# BY AGE & SEX (not really used; supplementary) #####
usoc$ageXsex <- as.factor(paste0(usoc$sex_dv, usoc$age.group))
means.byageXsex <- c()
for (group in c("Female10-12", "Female13-15", "Male10-12" , "Male13-15")) {
all.res <- c() # holds all results for one group (e.g. females)
for (year.x in sort(unique(usoc$int.year.sched)) ){
xs <- usoc[usoc$int.year.sched==year.x,] # isolate data for this YEAR
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
tmp <- usoc.svy$variables$ageXsex == group # limit survey data to group
usoc.svg.group <- subset(usoc.svy, tmp)
options(survey.lonely.psu = "adjust") # tells the survey package to handle singleton PSUs by scaling their variance estimates
group.glm <- lapply(dvlist, function(x) {
svyglm(eval(substitute(i ~ 1,
list(i = as.name(x)))),
design = usoc.svg.group)})
# extract cross-sectional results for all domains in this YEAR and for this GROUP
xs.res <- c()
for (domain in 1:6){
xs.res <- as.data.frame(rbind(xs.res,
cbind(year.x, group, dvlist[domain],
t(summary(group.glm[[domain]])[["coefficients"]][1:2]))))
}
names(xs.res) <- c("year.x", "group", "dvlist", "mean", "se")
all.res <- rbind(all.res, xs.res) # merge results from all years
all.res$mean <- as.numeric(all.res$mean)
all.res$se <- as.numeric(all.res$se)
} # finishes analysis for all years
means.byageXsex <- rbind(means.byageXsex, all.res)
} # finishes analyses for both groups
means.byageXsex$year.x <- as.numeric(means.byageXsex$year.x)
means.byageXsex$group_nice <- str_replace(means.byageXsex$group, "ale", "ales aged ")
# BY INCOME #####
means.byincome <- c()
props.byincome <-c()
for (income in c("high.inc", "mid.inc" , "low.inc")) {
all.res <- c()
for (year.x in sort(unique(usoc$int.year.sched)) ){
xs <- usoc[usoc$int.year.sched==year.x , ] # isolate data for one year
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
tmp <- usoc.svy$variables$income.groups == income # limit survey data to group
usoc.svg.group <- subset(usoc.svy, tmp)
options(survey.lonely.psu = "adjust") # tells the survey package to handle singleton PSUs by scaling their variance estimates
group.glm <- lapply(dvlist, function(x) {
svyglm(eval(substitute(i ~ 1,
list(i = as.name(x)))),
design = usoc.svg.group)})
# extract cross-sectional results for all domains in this YEAR and for this GROUP
xs.res <- c()
for (domain in 1:6){
xs.res <- as.data.frame(rbind(xs.res,
cbind(year.x, income, dvlist[domain],
t(summary(group.glm[[domain]])[["coefficients"]][1:2]))))
}
names(xs.res) <- c("year.x", "group", "dvlist", "mean", "se")
all.res <- rbind(all.res, xs.res)
# find weighted proportion of unhappy faces in this year (with 95% CI)
prop <- lapply(dvlist2, function(x) {
svyciprop(eval(substitute(i ~ 1 ,
list(i = as.name(x)))),
design = usoc.svg.group)})
y <- c()
for (domain in 1:6){
y <- as.data.frame(rbind(y,
data.frame( cbind(year.x,
income,
names(prop[[domain]]),
prop[[domain]],
confint(prop[[domain]]) ) )))
}
rownames(y) <- c()
names(y) <- c("year.x", "group", "domain", "prop", "ci.low", "ci.hi")
y$year.x <- as.numeric(y$year.x)
y$prop <- as.numeric(y$prop)
y$ci.low <- as.numeric(y$ci.low)
y$ci.hi <- as.numeric(y$ci.hi)
}
means.byincome <- rbind(means.byincome, all.res)
props.byincome <- rbind(means.byincome, all.res)
} # end of calculating means for each wave for each income group
means.byincome$year.x <- as.numeric(means.byincome$year.x)
means.byincome$mean <- as.numeric(means.byincome$mean)
means.byincome$se <- as.numeric(means.byincome$se)
#means.byincome$group <- factor(means.byincome$group, levels=c(3,2,1), labels=c("high.inc","mid.inc", "low.inc"))
# BY COUNTRY #####
means.bycountry <- c()
props.bycountry <- c()
for (country in c("England" , "Wales", "Scotland" ,"Northern Ireland")) {
all.group.means <- c() # holds all mean ratings for one group (e.g. females)
all.group.props <- c() # does same but for proportions
for (year.x in sort(unique(usoc$int.year.sched)) ){
if (year.x == 2024 & country=="Northern Ireland") { # n=0 in northern ireland in 2024
# do nothing
} else {
xs <- usoc[usoc$int.year.sched==year.x,] # isolate data for this YEAR
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
tmp <- usoc.svy$variables$country == country # limit survey data to group
usoc.svg.group <- subset(usoc.svy, tmp)
group.n <- nrow(usoc.svg.group$variables)
options(survey.lonely.psu = "adjust") # tells the survey package to handle singleton PSUs by scaling their variance estimates
group.glm <- lapply(dvlist, function(x) {
svyglm(eval(substitute(i ~ 1,
list(i = as.name(x)))),
design = usoc.svg.group)})
# extract cross-sectional results for all domains in this YEAR and for this GROUP
xs.res <- c()
for (domain in 1:6){
xs.res <- as.data.frame(rbind(xs.res,
cbind(year.x,
country,
dvlist[domain],
group.n,
t(summary(group.glm[[domain]])[["coefficients"]][1:2]))))
}
# next find weighted proportion of unhappy faces in this year, for this group (with 95% CI)
prop <- lapply(dvlist2, function(x) {
svyciprop(eval(substitute(i ~ 1 ,
list(i = as.name(x)))),
design = usoc.svg.group)})
y <- c()
for (domain in 1:6){
y <- as.data.frame(rbind(y,
data.frame( cbind(year.x,
country,
group.n,
names(prop[[domain]]),
prop[[domain]],
confint(prop[[domain]]) ) )))
}
rownames(y) <- c()
names(y) <- c("year.x","group", "sample.n", "domain", "prop", "ci.low", "ci.hi")
y$year.x <- as.numeric(y$year.x)
y$prop <- as.numeric(y$prop)
y$ci.low <- as.numeric(y$ci.low)
y$ci.hi <- as.numeric(y$ci.hi)
all.group.props <- rbind(all.group.props, y)
names(xs.res) <- c("year.x", "group", "dvlist","group.n", "mean", "se")
all.group.means <- rbind(all.group.means, xs.res) # merge results from all years
all.group.means$mean <- as.numeric(all.group.means$mean)
all.group.means$se <- as.numeric(all.group.means$se)
}} # finishes analysis of one year
means.bycountry <- rbind(means.bycountry, all.group.means)
props.bycountry <- rbind(props.bycountry, all.group.props)
} # finishes analyses for all groups
means.bycountry$year.x <- as.numeric(means.bycountry$year.x)
save(means.fullsample, prop.fullsample, means.byage, means.byage2, means.byage3, means.bysex, means.byageXsex, means.byincome,
props.byage, props.byage2, props.bysex, props.byincome, means.bycountry, props.bycountry,
file="/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/cross-sec_means.Rdata")
rm(list=c( "group.glm", "all.res", "usoc.svg.group", "xs", "xs.res", "usoc.svy", "x" ))# Plot distribution of satisfaction, changing over time, with weighted mean
dat <- c()
for (year in sort(unique(usoc$int.year.sched))){
d <- data.frame(year, table(usoc$h.life.c[usoc$int.year.sched==year]))
names(d) <- c("year", "response", "n")
d$tot.n <- sum(d$n) # note: this is the N who answered the Q
d$mode <- as.factor(ifelse(d$response == d$response[d$n==max(d$n)], 1, 0))
# get model-derived mean (pop weighted)
d$mean <- means.fullsample$mean[means.fullsample$domain=="h.life.c" & means.fullsample$year.x==year]
dat <- rbind(dat, d)
}
dat$prop <- dat$n/dat$tot.n
dat$response <- as.numeric(dat$response)-4
ggplot(data = dat)+
geom_col(aes(x=response, y=prop, fill = mode))+
geom_vline(aes(xintercept = mean), colour = "red")+
geom_label(aes(x=mean, y = .10, label= paste0("M = ",round(mean,2) )), size = 2, colour="red")+
scale_y_continuous(labels = scales::percent_format(accuracy = 5L), name="% of sample that year")+
scale_x_continuous(breaks=c(-3,0,3), labels=c(":-(", ":-|", ":-)"), name = "\nVisual response scale (coded -3 to +3)")+
labs(title = "The shifting distribution of life satisfaction (2009-2024)")+
scale_fill_manual(values = c("grey40", "black"))+
facet_wrap(~year)+
theme(plot.title = element_text(hjust = 0.5, face="bold"))
Calculate change statistics (linear trends)
It’s at this point that we go from weighted cross-sectional estimates, to trends. First, we simply place the cross-sectional points on a timeline and join them up graphically (linearly, using a LOESS line of best fit, etc.). Second, we run a GLM that assesses the slope of the linear change across all these means (and within 5-year windows). Use a model that acknowledges the inter-relatedness of data points, especially neighbouring points.
# TEST SIGNIFICANCE OF LINEAR TRENDS in the population mean
# save meta-regression output here:
trend.avgs <- c()
trend.full <- c()
trends.windows <- c()
trends.within.windows <- c() # not used
for (domain in c("h.school.c","h.schoolwork.c", "h.friends.c","h.fam.c", "h.appear.c","h.life.c" )) {
# limit data to cross-sectional means for THAT domain
dat <- means.fullsample[means.fullsample$domain==domain, ]
# create variance-covariance matrix for this data
# step 1: create diagnoal matrix of standard errors for this domain's mean
D <- diag(dat$SE, 16, 16)
# create variance-covar matrix by multiplying SEs (squared) by the correlation matrix R
# R describes how each timepoint is correlated with each other
V <- D %*% R %*% D
# run an intercept-only meta-regression: What is the pooled mean across all annual estimates?
rma.ave <- rma.mv(yi = mean,
V = V, # note: the standard errors for each annual mean are acknowledged in V
data = dat)
# save
trend.avgs <- rbind(trend.avgs, data.frame(
domain = domain,
beta = as.numeric(rma.ave[["beta"]]),
se = as.numeric(rma.ave[["se"]]),
ci.low = as.numeric(rma.ave[["ci.lb"]]),
ci.hi = as.numeric(rma.ave[["ci.ub"]])) )
# include time as moderator/covariate...
# full time period
rma.trend.full <- rma.mv(yi = mean,
mods = ~ year.n,
V = V,
data = dat)
# save slope estimates (change 2009-2024) for each domain
trend.full <- rbind(trend.full, data.frame(domain=domain,
beta=as.numeric(rma.trend.full[["beta"]][2]),
se=as.numeric(rma.trend.full[["se"]][2]),
ci.low=as.numeric(rma.trend.full[["ci.lb"]][2]),
ci.hi=as.numeric(rma.trend.full[["ci.ub"]][2]),
p=round(as.numeric(rma.trend.full[["pval"]][2]),4)))
# compare slopes within 5 year windows
dat$windows <- relevel(dat$windows, "2")
rma.trends.inwin <- rma.mv(yi = mean,
mods = ~ year.n*windows,
V = V,
data = dat)
# save slope estimates (change 2009-2024) for each domain (NOT USED)
trends.within.windows <- rbind(trends.within.windows,
data.frame(domain=domain,
effect = rownames(rma.trends.inwin[["beta"]])[5],
beta=as.numeric(rma.trends.inwin[["beta"]][5]),
se=as.numeric(rma.trends.inwin[["se"]][5]),
ci.low=as.numeric(rma.trends.inwin[["ci.lb"]][5]),
ci.hi=as.numeric(rma.trends.inwin[["ci.ub"]][5]),
p=round(as.numeric(rma.trends.inwin[["pval"]][5]),4)),
data.frame(domain=domain,
effect = rownames(rma.trends.inwin[["beta"]])[6],
beta=as.numeric(rma.trends.inwin[["beta"]][6]),
se=as.numeric(rma.trends.inwin[["se"]][6]),
ci.low=as.numeric(rma.trends.inwin[["ci.lb"]][6]),
ci.hi=as.numeric(rma.trends.inwin[["ci.ub"]][6]),
p=round(as.numeric(rma.trends.inwin[["pval"]][6]),4)))
for (time.window in 1:3){
dat.win <- dat[!is.na(dat$windows) & dat$windows==time.window & !is.na(dat$windows), ]
D <- diag(dat.win$SE)
# limit the covariance matrix R to the correct years
Rsub <- R[2:nrow(R), 2:ncol(R)] # remove 2009 (excluded from 5-y windows)
if (time.window==1){
Rsub <- Rsub[1:5, 1:5]
} else if (time.window==2){
Rsub <- Rsub[6:10, 6:10]
} else if (time.window==3){
Rsub <- Rsub[11:15, 11:15]
}
# create the variance-covariance matrix
V <- D %*% Rsub %*% D
this.windows.rma <- rma.mv(yi = mean,
V=V,
mods = ~ year.n,
data = dat.win)
# save slope estimates (with this 5-year window) for each domain
this.windows.res <- data.frame(
window = paste0("win", time.window),
domain=domain,
beta = as.numeric(this.windows.rma[["beta"]][2]),
se = as.numeric(this.windows.rma[["se"]][2]),
ci.low = as.numeric(this.windows.rma[["ci.lb"]][2]),
ci.hi = as.numeric(this.windows.rma[["ci.ub"]][2]),
p = round(as.numeric(this.windows.rma[["pval"]][2]),4))
trends.windows <- rbind(trends.windows, this.windows.res)
} # go to next 5-year window until all 3 are completed
} # go to next domain until all 6 are completed
# extract slopes for each 5-year segment (for a heatmap table)
slopes <- as.data.frame(cbind(trends.windows$domain[trends.windows$window=="win1"],
round(trends.windows$beta[trends.windows$window=="win1"], 3),
round(trends.windows$beta[trends.windows$window=="win2"], 3),
round(trends.windows$beta[trends.windows$window=="win3"], 3) ))
names(slopes) <- c("domain", "win1", "win2", "win3")
slopes$order <- c(3,4,2,1,5,6)
trend.full[,c(2:5)] <- round(trend.full[,c(2:5)],3)
# apply another meta-regression to the slopes within each 5-year window (NOT USED)
trend.diffs <- c()
for (domain in c("h.school.c","h.schoolwork.c", "h.friends.c","h.fam.c", "h.appear.c","h.life.c" )) {
# limit data to cross-sectional means for THAT domain
dat <- trends.windows[trends.windows$domain==domain, ]
# try another approach: simple z-test (USED)
z_win1win2 = dat$beta[1] - dat$beta[2] / sqrt(dat$se[1]^2 + dat$se[2]^2)
p_win1win2 = 2*pnorm(abs(z_win1win2), lower.tail = F)
z_win2win3 = dat$beta[2] - dat$beta[3] / sqrt(dat$se[2]^2 + dat$se[3]^2)
p_win2win3 = 2*pnorm(abs(z_win2win3), lower.tail = F)
# save slope estimates (with this 5-year window) for each domain
trend.diffs <- rbind(trend.diffs, data.frame(
domain=domain,
z_win1win2=z_win1win2, p_win1win2=round(p_win1win2,4),
z_win2win3=z_win2win3, p_win2win3=round(p_win2win3,4)))
}
print(trend.avgs)## domain beta se ci.low ci.hi
## 1 h.school.c 1.571920 0.01650794 1.539565 1.604275
## 2 h.schoolwork.c 1.365224 0.01399184 1.337801 1.392648
## 3 h.friends.c 2.301738 0.01056611 2.281029 2.322447
## 4 h.fam.c 2.374952 0.01105875 2.353277 2.396627
## 5 h.appear.c 1.263171 0.01586293 1.232080 1.294262
## 6 h.life.c 1.832274 0.01303964 1.806717 1.857831print(trend.full)## domain beta se ci.low ci.hi p
## 1 h.school.c -0.028 0.004 -0.035 -0.021 0.0000
## 2 h.schoolwork.c -0.019 0.003 -0.025 -0.013 0.0000
## 3 h.friends.c -0.025 0.002 -0.029 -0.020 0.0000
## 4 h.fam.c -0.006 0.002 -0.011 -0.002 0.0066
## 5 h.appear.c -0.018 0.004 -0.025 -0.011 0.0000
## 6 h.life.c -0.028 0.003 -0.033 -0.022 0.0000# what proportion of people reported being unhappy with life (sad faces) between 2009 and 2014?
mean(prop.fullsample$prop[prop.fullsample$domain=="unhappy.life" & (prop.fullsample$windows==1)], na.rm=T)## [1] 0.03932591#... between 2015 and 2019?
mean(prop.fullsample$prop[prop.fullsample$domain=="unhappy.life" & (prop.fullsample$windows==2)], na.rm=T)## [1] 0.05639085#... between 2020 and 2024?
mean(prop.fullsample$prop[prop.fullsample$domain=="unhappy.life" & (prop.fullsample$windows==3)], na.rm=T)## [1] 0.07472576# TEST SIGNIFICANCE OF LINEAR TRENDS in population variance
# save meta-regression output here:
var.avgs <- c()
var.16ytrend <- c()
var.5ytrends <- c()
for (domain in c("h.school.c","h.schoolwork.c", "h.friends.c","h.fam.c", "h.appear.c","h.life.c" )) {
# limit data to cross-sectional variance for THAT domain
dat <- means.fullsample[means.fullsample$domain==domain, ]
# create variance-covariance matrix for this data
# step 1: create diagnoal matrix of standard errors for this domain's mean
D <- diag(dat$popVar.SE, 16, 16)
# create variance-covar matrix by multiplying SEs (squared) by the correlation matrix R
# R describes how each timepoint is correlated with each other
V <- D %*% R %*% D
# run an intercept-only M.A: What is the pooled mean across all annual estimates?
rma.ave <- rma.mv(yi = popVar,
V = V, # note: the standard errors for each annual mean are acknowledged in V
data = dat)
# save
var.avgs <- rbind(var.avgs, data.frame(
domain = domain,
beta = as.numeric(rma.ave[["beta"]]),
se = as.numeric(rma.ave[["se"]]),
ci.low = as.numeric(rma.ave[["ci.lb"]]),
ci.hi = as.numeric(rma.ave[["ci.ub"]])) )
# include time as moderator...
# full time period
rma.trend.full <- rma.mv(yi = popVar,
mods = ~ year.n,
V = V,
data = dat)
# save slope estimates (change 2009-2024) for each domain
var.16ytrend <- rbind(var.16ytrend, data.frame(domain=domain,
beta=as.numeric(rma.trend.full[["beta"]][2]),
se=as.numeric(rma.trend.full[["se"]][2]),
ci.low=as.numeric(rma.trend.full[["ci.lb"]][2]),
ci.hi=as.numeric(rma.trend.full[["ci.ub"]][2]),
p=round(as.numeric(rma.trend.full[["pval"]][2]),4)))
# note the script below is pulling out the change in variance over time (within each 5y window)
# not the mean variance in a 5-year window
for (time.window in 1:3){
dat.win <- dat[!is.na(dat$windows) & dat$windows==time.window & !is.na(dat$windows), ]
D <- diag(dat.win$popVar.SE)
# limit the covariance matrix R to the correct years
Rsub <- R[2:nrow(R), 2:ncol(R)] # remove 2009 (excluded from 5-y windows)
if (time.window==1){
Rsub <- Rsub[1:5, 1:5]
} else if (time.window==2){
Rsub <- Rsub[6:10, 6:10]
} else if (time.window==3){
Rsub <- Rsub[11:15, 11:15]
}
# create the variance-covariance matrix
V <- D %*% Rsub %*% D
this.windows.rma <- rma.mv(yi = popVar,
V=V,
mods = ~ year.n,
data = dat.win)
# save variance estimates (with this 5-year window) for each domain
this.windows.res <- data.frame(
window = paste0("win", time.window),
domain=domain,
beta = as.numeric(this.windows.rma[["beta"]][2]),
se = as.numeric(this.windows.rma[["se"]][2]),
ci.low = as.numeric(this.windows.rma[["ci.lb"]][2]),
ci.hi = as.numeric(this.windows.rma[["ci.ub"]][2]),
p = round(as.numeric(this.windows.rma[["pval"]][2]),4))
var.5ytrends <- rbind(var.5ytrends, this.windows.res)
} # go to next 5-year window until all 3 are completed
} # go to next domain until all 6 are completed
# extract slopes for each 5-year segment (for a heatmap table)
slopes <- as.data.frame(cbind(var.5ytrends$domain[var.5ytrends$window=="win1"],
var.5ytrends$beta[var.5ytrends$window=="win1"],
var.5ytrends$beta[var.5ytrends$window=="win2"],
var.5ytrends$beta[var.5ytrends$window=="win3"]))
names(slopes) <- c("domain", "win1", "win2", "win3")
slopes$order <- c(3,4,2,1,5,6)
var.16ytrend[,c(2:5)] <- round(var.16ytrend[,c(2:5)],3)# TEST SIGNIFICANCE OF LINEAR TRENDS
# calculate SE from 95% CIs
prop.fullsample$SE <- (prop.fullsample$ci.hi - prop.fullsample$prop)/1.96
prop.fullsample$SE2 <- (prop.fullsample$prop - prop.fullsample$ci.low)/1.96
# NOTE: NOT SYMMETRIC SE's-- what's best way to get SE from a CI? Use larger of the two SEs to be conservative
# save meta-regression output here:
trend.avgs.unhappy <- c()
trend.full.unhappy <- c()
trends.windows.unhappy <- c()
dvlist2 <- c("unhappy.school", "unhappy.schoolwk", "unhappy.friends","unhappy.fam", "unhappy.appear", "unhappy.life" )
for (domain in dvlist2) {
dat <- prop.fullsample[prop.fullsample$domain==domain,]
# constuct var-covar matrix
D <- diag(dat$SE, 16, 16)
V <- D %*% R %*% D
# intercept-only M.A: What is the pooled proportion across all annual estimates?
rma.ave <- rma.mv(yi = prop,
V=V,
data = dat)
# save
trend.avgs.unhappy <- rbind(trend.avgs.unhappy, data.frame(
domain = domain,
beta = as.numeric(rma.ave[["beta"]]),
se = as.numeric(rma.ave[["se"]]),
ci.low = as.numeric(rma.ave[["ci.lb"]]),
ci.hi = as.numeric(rma.ave[["ci.ub"]])) )
# include time as moderator...
# full time period
rma.trend.full <- rma.mv(yi = prop,
V=V,
mods = ~ year.n,
data = dat)
# save slope estimates (change 2009-2024) for each domain
trend.full.unhappy <- rbind(trend.full.unhappy, data.frame(
domain=domain,
beta=as.numeric(rma.trend.full[["beta"]][2]),
se=as.numeric(rma.trend.full[["se"]][2]),
ci.low=as.numeric(rma.trend.full[["ci.lb"]][2]),
ci.hi=as.numeric(rma.trend.full[["ci.ub"]][2]),
p=round(as.numeric(rma.trend.full[["pval"]][2]),4)))
for (time.window in 1:3){
# restrict data again to this 5-year time window
dat.win <- dat[!is.na(dat$windows) & dat$windows==time.window, ]
D <- diag(dat.win$SE)
# limit the covariance matrix R to the correct years
Rsub <- R[2:nrow(R), 2:ncol(R)] # remove 2009 (excluded from 5-y windows)
if (time.window==1){
Rsub <- Rsub[1:5, 1:5]
} else if (time.window==2){
Rsub <- Rsub[6:10, 6:10]
} else if (time.window==3){
Rsub <- Rsub[11:15, 11:15]
}
# create the variance-covariance matrix
V <- D %*% Rsub %*% D
this.windows.rma <- rma.mv(yi = prop,
V=V,
mods = ~ year.n,
data = dat.win)
# save slope estimates (with this 5-year window) for each domain
this.windows.res <- data.frame(
window = paste0("win", time.window),
domain=domain,
beta = as.numeric(this.windows.rma[["beta"]][2]),
se = as.numeric(this.windows.rma[["se"]][2]),
ci.low = as.numeric(this.windows.rma[["ci.lb"]][2]),
ci.hi = as.numeric(this.windows.rma[["ci.ub"]][2]),
p = round(as.numeric(this.windows.rma[["pval"]][2]),4))
trends.windows.unhappy <- rbind(trends.windows.unhappy, this.windows.res)
} # go to next 5-year window until all 3 are completed
} # go to next domain until all 6 are completed
# extract slopes for each 5-year segment (for a heatmap table)
slopes.unhappy <- as.data.frame(cbind(trends.windows.unhappy$domain[trends.windows.unhappy$window=="win1"],
round(trends.windows.unhappy$beta[trends.windows.unhappy$window=="win1"]*100, 3),
round(trends.windows.unhappy$beta[trends.windows.unhappy$window=="win2"]*100, 3),
round(trends.windows.unhappy$beta[trends.windows.unhappy$window=="win3"]*100, 3)))
names(slopes.unhappy) <- c("domain", "win1", "win2", "win3")
slopes.unhappy$order <- c(3,4,2,1,5,6)
# express beta coefficients as % (which is the unit of the outcome)
trend.full.unhappy$beta.perc <- round(trend.full.unhappy$beta*100,2)
trend.full.unhappy$se.perc <- round(trend.full.unhappy$se*100,2)
trend.full.unhappy$ci.low.perc <- round(trend.full.unhappy$ci.low*100,2)
trend.full.unhappy$ci.hi.perc <- round(trend.full.unhappy$ci.hi*100,2)# run meta-regressions (custom function built in UndSocfunctions.r) to identify statistically significant differences exist between groups...
means.byage$group <- relevel(as.factor(means.byage$group), "13-15")
means.byincome$group <- relevel(as.factor(means.byincome$group), "mid.inc")
# make older group the ref so coefficients are mostly positive
# do slopes across the 15 years vary significantly between sexes, ages and hsd incomes? also shows group differences in means across the full time period.
means.bysex.15ytrend <- metareg.timeXgroup(means.bysex, "full")
means.byage.15ytrend <- metareg.timeXgroup(means.byage, "full")
means.byage3$group <- relevel(as.factor(means.byage3$group), "12-13")
means.byage3.15ytrend <- metareg.timeXgroup(means.byage3, "full")
means.byinc.15ytrend <- metareg.timeXgroup(means.byincome, "full")
#means.bysex.5y <- metareg.timeXgroup(means.bysex, "5years")
# try get this working to see when sex differences in trends began to emergePlots of trends
# plot means over time for full sample
col.palette = cols6
means.fullsample$`Happiness with...` <- factor(means.fullsample$domain,
levels=c( "h.life.c", "h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
# for (x in unique(means.fullsample$domain)){
# means.fullsample$lin.slope[means.fullsample$domain==x] <- round(trend.full$beta[trend.full$domain==x],2)
# }
# MEANS
# plot population-weighted trends (mean)
ggplot(data=means.fullsample,
aes(x=year.x, y=mean, group=`Happiness with...`, fill=`Happiness with...`))+
geom_point(aes(colour=`Happiness with...`))+
geom_line(aes(colour=`Happiness with...`), size=1.5)+
geom_ribbon(aes(ymin=mean-2*SE, ymax=mean+2*SE), alpha=0.2)+
#geom_text(aes(label=lin.slope, group=`Happiness with...`, x=2023, y=mean))+
labs(y="Mean happiness (-3 to +3)\n")+
scale_colour_manual(values=col.palette)+
scale_fill_manual(values=col.palette)+
scale_x_continuous(breaks=seq(2010,2024,2), minor_breaks=seq(2009,2023,2), limits = c(2009,2024), name = "Year")+
scale_y_continuous(breaks=c(1,2,3), minor_breaks = c(0.5, 1.5,2.5), limits = c(0.5, 2.5))+
theme_minimal()+
theme(axis.text.x = element_text(angle=40, hjust = 1, size=11),
axis.text.y = element_text(size=11),
axis.title.y = element_text(size=12))
# fit trend lines
#.. loess
ggplot(data=means.fullsample, aes(x=year.x, y=mean, group=`Happiness with...`, fill=`Happiness with...`))+
geom_point(aes(colour=`Happiness with...`))+
geom_smooth(aes(colour=`Happiness with...`), method="loess")+
labs(y="\nMean happiness (-3 to +3)")+
scale_colour_manual(values=col.palette)+
scale_fill_manual(values=col.palette)+
scale_y_continuous(breaks=c(1,2,3), minor_breaks = c(1.5,2.5), limits = c(0.8,2.5))+
theme_minimal()+
theme(axis.text.x = element_text(angle=30, hjust = 1))
# ...linear
ggplot(data=means.fullsample, aes(x=year.x, y=mean, group=`Happiness with...`, fill=`Happiness with...`))+
geom_point(aes(colour=`Happiness with...`))+
geom_smooth(aes(colour=`Happiness with...`), method="lm", formula = y~x)+
labs(y="Mean happiness (-3 to +3)")+
scale_colour_manual(values=col.palette)+
scale_fill_manual(values=col.palette)+
scale_y_continuous(breaks=c(1,2,3), minor_breaks = c(1.5,2.5), limits = c(0.8,2.5))+
theme_minimal()+
theme(axis.text.x = element_text(angle=30, hjust = 1))
# ... linear within windows
means.fullsample$domainXwindow <- paste0(means.fullsample$`Happiness with...`, means.fullsample$windows)
means.fullsample$domainXwindow <- factor(means.fullsample$domainXwindow,
levels=c(paste0("Family",seq(1,3)),
paste0("Friends",seq(1,3)),
paste0("School",seq(1,3)),
paste0("Schoolwork",seq(1,3)),
paste0("Appearance",seq(1,3)),
paste0("Life as whole", seq(1,3))))
ggplot(data=means.fullsample, aes(x=year.x, y=mean, group=domainXwindow, fill=domainXwindow))+
geom_point(aes(colour=domainXwindow))+
geom_smooth(aes(colour=domainXwindow), method="lm", formula = y~x)+
labs(y="\nMean happiness (-3 to +3)", x="5-year windows (2009-2013; 2014-2018; 2019-2023)")+
scale_colour_manual(values=cols18)+
scale_fill_manual(values=cols18)+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
scale_y_continuous(breaks=c(1,2,3), minor_breaks = c(1.5,2.5), limits = c(0.5,2.5))+
theme_minimal()
Precision
It seems my mean estimates are becoming less precise over time. This is
likely driven by reductions in the sample size. It could also be driven
by reduced efficiency of weights or differences in subgroup sizes. But
we are assured by the UKHLS team that this is not the case. This is the
uncertainty around the mean estimate.
Spread
We don’t know from the error around the mean though whether the spread
or population variance in satisfaction is widening over time. That would
suggest widening differences between subgroups in the population.
# plot population variance over time
col.palette = cols6
means.fullsample$`Satisfaction with...` <- factor(means.fullsample$domain,
levels=c( "h.life.c", "h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
# for (x in unique(means.fullsample$domain)){
# means.fullsample$lin.slope[means.fullsample$domain==x] <- round(trend.full$beta[trend.full$domain==x],2)
# }
ggplot(data=means.fullsample,
aes(x=year.x, y=popVar, group=`Satisfaction with...`, fill=`Satisfaction with...`))+
geom_point(aes(colour=`Satisfaction with...`))+
geom_ribbon(aes(ymin=popVar-2*popVar.SE, ymax=popVar+2*popVar.SE), alpha=0.2)+
geom_line(aes(colour=`Satisfaction with...`), size=1.5)+
labs(y="Variance across the population")+
scale_colour_manual(values=col.palette)+
scale_fill_manual(values=col.palette)+
scale_x_continuous(breaks=seq(2010,2024,2), minor_breaks=seq(2009,2023,2), limits = c(2009,2024), name = "Year")+
#scale_y_continuous(breaks=c(1,2,3), minor_breaks = c(0.5, 1.5,2.5), limits = c(0.5, 2.5))+
theme_minimal()+
theme(axis.text.x = element_text(angle=40, hjust = 1, size=11),
axis.text.y = element_text(size=11),
axis.title.y = element_text(size=12))
Note that the y-axis units here are the likert satisfaction scale -3 to
+3
# trends in proportion unhappy
col.palette = cols6
prop.fullsample$`Unhappiness with...` <- factor(prop.fullsample$domain,
levels=c( "unhappy.life", "unhappy.fam", "unhappy.friends",
"unhappy.school" , "unhappy.schoolwk" ,
"unhappy.appear" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
ggplot(data=prop.fullsample,
aes(x=year.x, y=prop, group=`Unhappiness with...`, fill=`Unhappiness with...`))+
geom_point(aes(colour=`Unhappiness with...`))+
geom_line(aes(colour=`Unhappiness with...`), size=1.5)+
geom_ribbon(aes(ymin=ci.low, ymax=ci.hi), alpha=0.2)+
#geom_text(aes(label=lin.slope, group=`Unhappiness with...`, x=2023, y=mean))+
labs(y="Mean happiness (-3 to +3)\n")+
scale_colour_manual(values=col.palette)+
scale_fill_manual(values=col.palette)+
scale_x_continuous(breaks=seq(2010,2024,2), minor_breaks=seq(2009,2023,2), limits = c(2009,2024), name = "Year")+
scale_y_continuous(labels = scales::percent_format(accuracy = 5L), minor_breaks=NULL, breaks= c(0,.05, .10, .15, .20), limits = c(0, .23), name="Proportion unhappy")+
#scale_y_continuous(breaks=c(1,2,3), minor_breaks = c(1.5,2.5), limits = c(0.5, 2.5))+
theme_minimal()+
theme(axis.text.x = element_text(angle=40, hjust = 1, size=11),
axis.text.y = element_text(size=11),
axis.title.y = element_text(size=12))
# plot as bars, in rows
ggplot(data=prop.fullsample, aes(x=year.x, y=prop, fill=`Unhappiness with...`))+
geom_col(alpha=0.7)+
geom_errorbar(aes(ymin=ci.low, ymax=ci.hi,colour=`Unhappiness with...`), width=0)+
labs(y="Mean happiness (-3 to +3)\n")+
scale_colour_manual(values=col.palette)+
scale_fill_manual(values=col.palette)+
scale_x_continuous(breaks=seq(2010,2024,2), minor_breaks=seq(2009,2023,2), limits = c(2008.5,2024.5), name = "Year")+
scale_y_continuous(labels = scales::percent_format(accuracy = 5L), minor_breaks=NULL, breaks= c(0,.05, .10, .15, .20), limits = c(0, .23), name="Proportion unhappy")+
facet_grid(rows=vars(`Unhappiness with...`))+
theme_minimal()+
theme(axis.text.x = element_text(angle=40, hjust = 1, size=9),
axis.text.y = element_text(size=9),
axis.title.y = element_text(size=12),
legend.position = "none")
# plot as bars, in grid
ggplot(data=prop.fullsample, aes(x=year.x, y=prop, fill=`Unhappiness with...`))+
geom_col(alpha=0.7)+
geom_errorbar(aes(ymin=ci.low, ymax=ci.hi,colour=`Unhappiness with...`), width=0)+
#geom_smooth(method="lm", formula=y~x, colour="black", fill="black", size=0.7)+
labs(y="Mean happiness (-3 to +3)\n")+
scale_colour_manual(values=col.palette)+
scale_fill_manual(values=col.palette)+
scale_x_continuous(breaks=seq(2010,2024,2), minor_breaks=seq(2009,2023,2), limits = c(2008.5,2024.5), name = "Year")+
scale_y_continuous(labels = scales::percent_format(accuracy = 5L), minor_breaks=NULL, breaks= c(0,.05, .10, .15, .20), limits = c(0, .23), name="Proportion unhappy")+
facet_wrap(~`Unhappiness with...`)+
theme_minimal()+
theme(axis.text.x = element_text(angle=40, hjust = 1, size=10),
axis.text.y = element_text(size=10),
axis.title.y = element_text(size=14),
axis.title.x = element_text(size=13),
strip.text.x = element_text(size=13),
legend.position = "none")
# plot trends for each demographic group
# prepare dataframes for plotting
means.bysex$`Happiness with...` <- factor(means.bysex$dvlist,
levels=c( "h.life.c" ,"h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
means.bysex$domainXsex <- as.factor(paste(means.bysex$dvlist, means.bysex$sex, sep="_"))
means.byage$`Happiness with...` <- factor(means.byage$dvlist,
levels=c( "h.life.c" ,"h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
means.byage$age <- paste(means.byage$group, "years")
means.byincome$`Happiness with...` <- factor(means.byincome$dvlist,
levels=c( "h.life.c" ,"h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
means.byincome$group <- factor(means.byincome$group, levels=c("high.inc", "mid.inc", "low.inc" ),
labels=c("High income", "Middle income", "Low income"))
LIMS = c(0.4, 2.7)
gsex <- ggplot(data=means.bysex, aes(x=year.x, y=mean, group=group))+
geom_line(aes(colour=group), size=1)+
geom_ribbon(aes(ymin=mean-2*se, ymax=mean+2*se, fill=group),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
scale_y_continuous(breaks=c(0,1,2), minor_breaks = c(0.5, 1.5,2.5),
limits = LIMS, name = " ")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
scale_shape_manual(values=c(3,15))+
scale_colour_manual(values=c("darkorange2", "turquoise3"))+
scale_fill_manual(values=c( "darkorange2", "turquoise3"))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=6))
gage <- ggplot(data=means.byage, aes(x=year.x, y=mean, group=group))+
geom_line(aes(colour=group), size=1)+
geom_ribbon(aes(ymin=mean-2*se, ymax=mean+2*se, fill=group),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
scale_y_continuous(breaks=c(0,1,2), minor_breaks = c(0.5,1.5,2.5),
limits = LIMS, name = " ")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
scale_shape_manual(values=c(3,15))+
scale_colour_manual(values=c( "darkseagreen4","olivedrab3"))+
scale_fill_manual(values=c( "darkseagreen4", "olivedrab3"))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=6))
gincome <- ggplot(data=means.byincome, aes(x=year.x, y=mean, group=group))+
geom_line(aes(colour=group), size=0.7)+
geom_ribbon(aes(ymin=mean-2*se, ymax=mean+2*se, fill=group),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
scale_y_continuous(breaks=c(0,1,2), minor_breaks = c(0.5,1.5,2.5),
limits = LIMS, name = " ")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
scale_shape_manual(values=c(3,15))+
scale_colour_manual(values=c( "darkslateblue","plum4","thistle" ))+
scale_fill_manual(values=c( "darkslateblue","plum4","thistle" ))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=6))
# "salmon4","tomato3","goldenrod1"
ggarrange(gsex, gage, gincome, nrow=3, ncol=1, legend = "none")
means.byage3$`Happiness with...` <- factor(means.byage3$dvlist,
levels=c( "h.life.c" ,"h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
ggplot(data=means.byage3, aes(x=year.x, y=mean, group=group))+
geom_line(aes(colour=group), size=1)+
geom_ribbon(aes(ymin=mean-2*se, ymax=mean+2*se, fill=group),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
scale_y_continuous(breaks=c(0,1,2,3), minor_breaks = c(0.5,1.5,2.5),
limits = c(-0.1, 2.9), name = "Mean satisfaction")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
coord_cartesian(ylim=c(0.5,2.9))+
scale_colour_manual(values=c("olivedrab3", "darkseagreen3","darkgreen" ))+
scale_fill_manual(values=c("olivedrab3", "darkseagreen3","darkgreen"))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=7))
means.bycountry$`Happiness with...` <- factor(means.bycountry$dvlist,
levels=c( "h.life.c" ,"h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
means.bycountry$group <- factor(means.bycountry$group, levels=c("England", "Scotland", "Wales", "Northern Ireland" ))
ggplot(data=means.bycountry, aes(x=year.x, y=mean, group=group))+
#geom_line(aes(colour=group), size=1)+
geom_smooth(aes(colour=group), formula=y~x, se=F)+
#geom_ribbon(aes(ymin=mean-2*se, ymax=mean+2*se, fill=group),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
# scale_y_continuous(breaks=c(0,1,2,3), minor_breaks = c(0.5,1.5,2.5),
# limits = c(-0.1, 2.9), name = "Mean satisfaction")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2023), name = "")+
labs(y="mean satisfaction (-3 to +3)")+
coord_cartesian(ylim=c(0.5,2.9))+
scale_colour_manual(values=c("tomato2", "dodgerblue4", "seagreen3", "orange"))+
scale_fill_manual(values=c( "tomato2", "dodgerblue4", "seagreen3", "orange"))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=7))
# sample sizes are too small-- trends are v noisy. Which is why I've shown them smoothed here without error bars.
# you could calculate 2-5 year averages and re-run plots means.byageXsex$`Happiness with...` <- factor(means.byageXsex$dvlist,
levels=c( "h.life.c" ,"h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Life as whole", "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
means.byageXsex$ci.low <- means.byageXsex$mean - means.byageXsex$se*1.96
means.byageXsex$ci.hi <- means.byageXsex$mean + means.byageXsex$se*1.96
ggplot(data=means.byageXsex, aes(x=year.x, y=mean, group=group_nice))+
geom_line(aes(colour=group_nice), size=1)+
geom_ribbon(aes(ymin=mean-2*se, ymax=mean+2*se, fill=group_nice),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
scale_y_continuous(breaks=c(0,1,2,3), minor_breaks = c(0.5,1.5,2.5),
limits = c(-0.1, 2.9), name = "Mean satisfaction")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
coord_cartesian(ylim=c(0.2,2.7))+
scale_colour_manual(values=c("tomato", "tomato4","cyan3","aquamarine4"))+
scale_fill_manual(values=c( "tomato", "tomato4","cyan3","aquamarine4"))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=7),
legend.title=element_blank())
# plot appearance trends only
ggplot(data=means.byageXsex[means.byageXsex$dvlist=="h.appear.c",], aes(x=year.x, y=mean, group=group_nice))+
geom_line(aes(colour=group_nice), size=1)+
geom_ribbon(aes(ymin=mean-2*se, ymax=mean+2*se, fill=group_nice),alpha=0.2)+
#facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
scale_y_continuous(breaks=c(0,1,2,3), minor_breaks = c(0.5,1.5,2.5),
limits = c(-0.1, 2.9), name = "Mean satisfaction")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
coord_cartesian(ylim=c(0.2,2.7))+
scale_colour_manual(values=c("tomato", "tomato4","cyan3","aquamarine4"))+
scale_fill_manual(values=c( "tomato", "tomato4","cyan3","aquamarine4"))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=12),
legend.title=element_blank())
2. Relative influence of sub-domains on general life satisfaction
Find inter-correlations between satisfaction in each aspect of life
# takes ~15 seconds
# create data frames to store B coef and SEs of domain onto general life satisfaction, for each wave of data
coefs.wgts <- c() # adjusted for all other domains of life
coefs.wgts.unadj <- c() # unadjusted
for (year.x in sort(unique(usoc$int.year.sched)) ){
#for (wave.x in c("a", "c", "e", "g", "i", "k", "m")){
# extract cross-sec data for wave x
xs <- usoc[usoc$int.year.sched==year.x,]
# exclude participants with missing data on any life satis var
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
# set svy design for this wave
usoc.svy <- svydesign(id=~psu + hidp, # hidp accounts for dependence among individuals within the same household by adjusting standard errors (via robust sandwich estimators)
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
options(survey.lonely.psu = "adjust")
# PREDICT GENERAL LIFE SATISFACTIONF FROM EACH DOMAIN (adjusting for all other domains)
svyglm.multi <- svyglm(formula = h.life.c ~ h.school.c + h.schoolwork.c + h.friends.c + h.fam.c + h.appear.c,
design = usoc.svy)
x= as.data.frame(summary(svyglm.multi)[["coefficients"]])
coefs.wgts <- rbind(coefs.wgts, c(as.numeric(coef(svyglm.multi)), # betas
as.numeric(x[,2])) ) # SEs
# PREDICT GENERAL LIFE SATISFACTIONF FROM EACH DOMAIN (bivariate-ly)
ivlist <- c("h.school.c", "h.schoolwork.c", "h.friends.c", "h.fam.c", "h.appear.c")
# run GLM substituting each area of life as only predictor
svyglm.uni <- lapply(ivlist, function(x) {
svyglm(eval(substitute(h.life.c ~ i ,
list(i = as.name(x)))),
design = usoc.svy)})
x <- c()
for (domain in 1:5){
x <- as.data.frame(rbind(x,
data.frame(year.x=year.x,
domain=ivlist[domain],
Beta = summary(svyglm.uni[[domain]])[["coefficients"]][2,1], # mean of the 1st (non-intercept) term
SE = summary(svyglm.uni[[domain]])[["coefficients"]][2,2], # this is the sampling error of the mean.
N = length(svyglm.uni[[1]][["effects"]])
))) }
coefs.wgts.unadj <- rbind(coefs.wgts.unadj, x)
}
coefs.wgts <- cbind(sort(unique(usoc$int.year.sched)), as.data.frame(round(coefs.wgts,4)))
names(coefs.wgts) <- c("year", "B.intercept", "B.school", "B.schoolwork", "B.friends", "B.family", "B.appearance",
"SE.intercept", "SE.school", "SE.schoolwork", "SE.friends", "SE.family", "SE.appearance")
coefs.wgts$year.n <- coefs.wgts$year-2009 # centers at 2009
coefs.wgts$window <- c(rep(1,5), rep(2,5), rep(3,6)) # note extra year in window3!
# re-structure to long-format data
long.wgts <- coefs.wgts %>%
dplyr::select(year, B.school, B.schoolwork, B.friends, B.family, B.appearance)%>%
gather(key = "domain", value="Beta", -year)
long.se <- coefs.wgts %>%
dplyr::select(year, SE.school, SE.schoolwork, SE.friends, SE.family, SE.appearance)%>%
gather(key = "domain", value="SE", -year)
wgts.long <- cbind(long.wgts, long.se[,"SE"])
names(wgts.long)[4] <- "SE"
wgts.long$`Happiness with...` <- factor(wgts.long$domain,
levels=c( "B.family" , "B.friends",
"B.school", "B.schoolwork",
"B.appearance" ),
labels=c( "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
wgts.long$year.n <- wgts.long$year-2009
wgts.long$windows <- cut(wgts.long$year, breaks=c(2009, 2014, 2019, 2025), labels=c(1, 2, 3))
# Plot change over time in asssociations
coefs.wgts.unadj$`Satisfaction with...` <- factor(coefs.wgts.unadj$domain,
levels=c( "h.fam.c" , "h.friends.c",
"h.school.c", "h.schoolwork.c",
"h.appear.c" ),
labels=c( "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
# bar chart UNADJUSTED
ggplot(data=coefs.wgts.unadj, aes(x=year.x, y=Beta))+
geom_col(aes(fill=`Satisfaction with...`), alpha=0.7)+
geom_errorbar(aes(ymin=Beta-SE, ymax=Beta+SE, colour=`Satisfaction with...`), width=0)+
geom_smooth(aes(x=year.x, y=Beta), method="loess", colour="grey60", se=F)+
geom_smooth(aes(x=year.x, y=Beta), method="lm", formula = y~x, colour="black", se=T)+
#geom_smooth(aes(x=year.x, y=Beta), method="lm", formula = y~poly(x,2), colour="black", se=F)+
labs(y="Unadjusted association with life satisfaction (B ± SE)\n", x="\nSurvey year")+
facet_grid(~`Satisfaction with...`)+
scale_x_continuous(breaks= seq(2010,2024,2), limits = c(2008,2025), name = "")+
scale_y_continuous(limits=c(0, 1))+
scale_fill_manual(values=cols5)+
scale_colour_manual(values=cols5)+
theme_minimal()+
theme(axis.text.x = element_text(size=8, angle=60, hjust=1, vjust = 1),
#panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())
# bar chart ADJUSTED
ggplot(data=wgts.long, aes(x=year, y=Beta))+
geom_col(aes(fill=`Happiness with...`), alpha=0.7)+
geom_errorbar(aes(ymin=Beta-SE, ymax=Beta+SE, colour=`Happiness with...`), width=0)+
geom_smooth(aes(x=year, y=Beta), method="loess", colour="grey60", se=F)+
geom_smooth(aes(x=year, y=Beta), method="lm", formula = y~x, colour="black", se=T)+
#geom_smooth(aes(x=year, y=Beta), method="lm", formula = y~poly(x,2), colour="black", se=F)+
labs(y="Adjusted association with life satisfaction (B ± SE)\n", x="\nSurvey year")+
facet_grid(~`Happiness with...`)+
scale_x_continuous(breaks= seq(2010,2024,2), limits = c(2008,2025), name = "")+
scale_y_continuous(limits=c(0, 0.5))+
scale_fill_manual(values=cols5)+
scale_colour_manual(values=cols5)+
theme_minimal()+
theme(axis.text.x = element_text(size=8, angle=60, hjust=1, vjust = 1),
#panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())
## WITHOUT COVID
# bar chart UNADJUSTED
ggplot(data=coefs.wgts.unadj[coefs.wgts.unadj$year.x!=2020 & coefs.wgts.unadj$year.x!=2021,], aes(x=year.x, y=Beta))+
geom_col(aes(fill=`Satisfaction with...`), alpha=0.7)+
geom_errorbar(aes(ymin=Beta-SE, ymax=Beta+SE, colour=`Satisfaction with...`), width=0)+
geom_smooth(aes(x=year.x, y=Beta), method="loess", colour="grey60", se=F)+
geom_smooth(aes(x=year.x, y=Beta), method="lm", formula = y~x, colour="black", se=T)+
#geom_smooth(aes(x=year.x, y=Beta), method="lm", formula = y~poly(x,2), colour="black", se=F)+
labs(y="UNADJUSTED association with life satisfaction (B ± SE)\n", x="\nSurvey year")+
facet_grid(~`Satisfaction with...`)+
scale_x_continuous(breaks= seq(2010,2024,2), limits = c(2008,2025), name = "")+
scale_y_continuous(limits=c(0, 1))+
scale_fill_manual(values=cols5)+
scale_colour_manual(values=cols5)+
theme_minimal()+
theme(axis.text.x = element_text(size=8, angle=60, hjust=1, vjust = 1),
#panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())
# bar chart ADJUSTED
ggplot(data=wgts.long[wgts.long$year!=2020 & wgts.long$year!=2021,], aes(x=year, y=Beta))+
geom_col(aes(fill=`Happiness with...`), alpha=0.7)+
geom_errorbar(aes(ymin=Beta-SE, ymax=Beta+SE, colour=`Happiness with...`), width=0)+
geom_smooth(aes(x=year, y=Beta), method="loess", colour="grey60", se=F)+
geom_smooth(aes(x=year, y=Beta), method="lm", formula = y~x, colour="black", se=T)+
#geom_smooth(aes(x=year, y=Beta), method="lm", formula = y~poly(x,2), colour="black", se=F)+
labs(y="Adjusted association with life satisfaction (B ± SE)\n", x="\nSurvey year")+
facet_grid(~`Happiness with...`)+
scale_x_continuous(breaks= seq(2010,2024,2), limits = c(2008,2025), name = "")+
scale_y_continuous(limits=c(0, 0.5))+
scale_fill_manual(values=cols5)+
scale_colour_manual(values=cols5)+
theme_minimal()+
theme(axis.text.x = element_text(size=8, angle=60, hjust=1, vjust = 1),
#panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())
Note the y-axis represents the association between each domain and “general life satisfaction”, after controlling for the other domains.
- Notable rise in the relative influence of “happiness with
appearance” on general life satisfcation, over full time period
- Rise in relative influence of “happiness with family”
- Smaller changes in friends, school and schoolwork
# save meta-regression output here:
wgt.avgs <- c()
wgt.trends <- c()
wgt.avg.5ywin <- c() # this one's different to the previous "change *within* 5 year windows". This time we're looking at 5-year *average* i.e. intercept
var.cov.matrix <- c()
var.cov.guide <- c()
for (domain in levels(as.factor(wgts.long$domain))) {
dat <- wgts.long[wgts.long$domain==domain, ]
D <- diag(dat$SE, 16, 16)
V <- D %*% R %*% D
# store these variance-covariance matrices for subgroup analyses later
var.cov.matrix <- abind(var.cov.matrix, V)
var.cov.guide <- cbind(var.cov.guide, domain)
# intercept-only M.A: What is the pooled mean across all my waves?
rma.ave <- rma.mv(yi = Beta,
V = V,
data=dat)
# save
wgt.avgs <- rbind(wgt.avgs, data.frame(
domain = domain,
beta = as.numeric(rma.ave[["beta"]]),
se = as.numeric(rma.ave[["se"]]),
ci.low = as.numeric(rma.ave[["ci.lb"]]),
ci.hi = as.numeric(rma.ave[["ci.ub"]])) )
# include time as moderator...
# full time period
rma.trend.full <- rma.mv(yi = Beta,
V = V,
mods = ~ year.n,
data = dat)
# save
wgt.trends <- rbind(wgt.trends, data.frame(
domain=domain,
beta=as.numeric(rma.trend.full[["beta"]][2]),
se=as.numeric(rma.trend.full[["se"]][2]),
ci.low=as.numeric(rma.trend.full[["ci.lb"]][2]),
ci.hi=as.numeric(rma.trend.full[["ci.ub"]][2]),
p=round(as.numeric(rma.trend.full[["pval"]][2]),4)))
# Note: There's so much year-to-year fluctuation with this data, I don't think we should be looking at within-5 year changes. But we could estimate *average* associations within each 5-year window, to average out some of those fluctuations...
for (win in 1:3){
dat.win <- wgts.long[wgts.long$domain==domain & wgts.long$windows==win & !is.na(wgts.long$windows), ]
D <- diag(dat.win$SE)
# limit the covariance matrix R to the correct years
Rsub <- R[2:nrow(R), 2:ncol(R)] # remove 2009 (excluded from 5-y windows)
if (win==1){
Rsub <- Rsub[1:5, 1:5]
} else if (win==2){
Rsub <- Rsub[6:10, 6:10]
} else if (win==3){
Rsub <- Rsub[11:15, 11:15]
}
# create the variance-covariance matrix
V <- D %*% Rsub %*% D
rma.ave.win <- rma.mv(yi = Beta, V = V, data=dat.win)
wgt.avg.5ywin <- rbind(wgt.avg.5ywin, data.frame(
domain = domain,
window = win,
beta = as.numeric(rma.ave.win[["beta"]]),
se = as.numeric(rma.ave.win[["se"]]),
ci.low = as.numeric(rma.ave.win[["ci.lb"]]),
ci.hi = as.numeric(rma.ave.win[["ci.ub"]])) )
}
}
wgt.avgs[,c(2:5)] <- round(wgt.avgs[,c(2:5)], 2)
wgt.trends[,c(2:5)] <- round(wgt.trends[,c(2:5)], 3)
print(wgt.avgs)## domain beta se ci.low ci.hi
## 1 B.appearance 0.25 0.01 0.24 0.27
## 2 B.family 0.33 0.01 0.31 0.36
## 3 B.friends 0.15 0.01 0.12 0.18
## 4 B.school 0.10 0.01 0.08 0.12
## 5 B.schoolwork 0.12 0.01 0.10 0.13print(wgt.trends)## domain beta se ci.low ci.hi p
## 1 B.appearance 0.009 0.002 0.005 0.012 0.0000
## 2 B.family 0.003 0.003 -0.003 0.008 0.3215
## 3 B.friends -0.004 0.003 -0.009 0.001 0.1453
## 4 B.school 0.000 0.002 -0.004 0.003 0.8751
## 5 B.schoolwork -0.001 0.002 -0.004 0.003 0.7309library(abind)
# save meta-regression output here:
wgt.avgs <- c()
wgt.trends <- c()
var.cov.matrix <- c()
var.cov.guide <- c()
for (domain in levels(as.factor(coefs.wgts.unadj$domain))) {
dat <- coefs.wgts.unadj[coefs.wgts.unadj$domain==domain, ]
D <- diag(dat$SE, 16, 16)
V <- D %*% R %*% D
# store these variance-covariance matrices for subgroup analyses later
var.cov.matrix <- abind(var.cov.matrix, V)
var.cov.guide <- cbind(var.cov.guide, domain)
# intercept-only M.A: What is the pooled mean across all my waves?
rma.ave <- rma.mv(yi = Beta,
V = V,
data=dat)
# save
wgt.avgs <- rbind(wgt.avgs, data.frame(
domain = domain,
beta = as.numeric(rma.ave[["beta"]]),
se = as.numeric(rma.ave[["se"]]),
ci.low = as.numeric(rma.ave[["ci.lb"]]),
ci.hi = as.numeric(rma.ave[["ci.ub"]])) )
# include time as moderator...
# full time period
rma.trend.full <- rma.mv(yi = Beta,
V = V,
mods = ~ year.x,
data = dat)
# save
wgt.trends <- rbind(wgt.trends, data.frame(
domain=domain,
beta=as.numeric(rma.trend.full[["beta"]][2]),
se=as.numeric(rma.trend.full[["se"]][2]),
ci.low=as.numeric(rma.trend.full[["ci.lb"]][2]),
ci.hi=as.numeric(rma.trend.full[["ci.ub"]][2]),
p=round(as.numeric(rma.trend.full[["pval"]][2]),4)))
}
wgt.avgs[,c(2:5)] <- round(wgt.avgs[,c(2:5)], 2)
wgt.trends[,c(2:5)] <- round(wgt.trends[,c(2:5)], 3)– by demographic group
# step 1. run the adjusted model (all subdomains predicting gen life happiness) within each subgroup
# TAKES ~1-2 MINS
wgts.bygroup <- c()
for (group in c("Female", "Male", "10-12", "13-15", "high.inc", "low.inc" , "mid.inc")){
for (year.x in sort(unique(usoc$int.year.sched)) ){
# limit dataset to a specific year
xs <- usoc[usoc$int.year.sched==year.x ,]
xs <- xs[!is.na(xs$h.school.c) & !is.na(xs$h.schoolwork.c) & !is.na(xs$h.friends.c) & !is.na(xs$h.fam.c) & !is.na(xs$h.appear.c) & !is.na(xs$h.life.c), ]
usoc.svy <- svydesign(id=~psu + hidp,
strata=~strata,
weights=~wgt.rescaled,
data=xs,
nest=TRUE)
options(survey.lonely.psu = "adjust")
# limit cross-sectional survey data to group
if (group=="Female"|group=="Male"){
tmp <- usoc.svy$variables$sex_dv == group
} else if (group=="10-12"|group=="13-15") {
tmp <- usoc.svy$variables$age.group == group
} else if (str_detect(group, "inc")) {
tmp <- usoc.svy$variables$income.groups == group
}
usoc.svg.group <- subset(usoc.svy, tmp)
options(survey.lonely.psu = "adjust")
# find relative influence of each aspect of life, for THIS group, in this year.
group.glm <- svyglm(formula = h.life.c ~ h.school.c + h.schoolwork.c
+ h.friends.c + h.fam.c + h.appear.c,
design = usoc.svg.group)
coef <- as.data.frame(summary(group.glm)[["coefficients"]])
wgts.bygroup <- rbind(wgts.bygroup, data.frame(
year = year.x,
group = group,
n = length(group.glm$effects),
school.B = coef$Estimate[rownames(coef)=="h.school.c"],
school.SE = coef$`Std. Error`[rownames(coef)=="h.school.c"],
schoolwork.B = coef$Estimate[rownames(coef)=="h.schoolwork.c"],
schoolwork.SE = coef$`Std. Error`[rownames(coef)=="h.schoolwork.c"],
fam.B = coef$Estimate[rownames(coef)=="h.fam.c"],
fam.SE = coef$`Std. Error`[rownames(coef)=="h.fam.c"],
friend.B = coef$Estimate[rownames(coef)=="h.friends.c"],
friend.SE = coef$`Std. Error`[rownames(coef)=="h.friends.c"],
appear.B = coef$Estimate[rownames(coef)=="h.appear.c"],
appear.SE = coef$`Std. Error`[rownames(coef)=="h.appear.c"]))
}
}
# re-structure to long-format data
long.wgts <- wgts.bygroup %>%
dplyr::select(year, group, n, school.B, schoolwork.B, friend.B, fam.B, appear.B)%>%
gather(key = "domain", value="Beta", -c(year, group, n))
long.se <- wgts.bygroup %>%
dplyr::select(year, group, n, school.SE, schoolwork.SE, friend.SE, fam.SE, appear.SE)%>%
gather(key = "domain", value="SE", -c(year, group, n))
group.wgts.long <- cbind(long.wgts, long.se[,"SE"])
names(group.wgts.long)[6] <- "SE"
group.wgts.long$`Happiness with...` <- factor(group.wgts.long$domain,
levels=c("fam.B","friend.B", "school.B", "schoolwork.B", "appear.B" ),
labels=c( "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
group.wgts.long$year.n <- group.wgts.long$year-2009
# step 2
# input group-specific estimates into a meta-analysis to estimate
#(i) the average association estimates for that group and
#(ii) see the association has changed over time for either group
group.B.ave <- c()
group.B.5yave <- c()
group.B.change <- c()
for (group in levels(as.factor(group.wgts.long$group))) {
for (domain in levels(as.factor(group.wgts.long$domain))) {
dat <- group.wgts.long[group.wgts.long$domain==domain & group.wgts.long$group==group, ]
D <- diag(dat$SE, 16, 16)
V <- D %*% R %*% D
# intercept-only M.A: What is the pooled estimate across all years?
rma.ave <- rma.mv(yi = Beta,
V = V,
data=dat)
# save
group.B.ave <- rbind(group.B.ave, data.frame(
domain = domain,
group = group,
beta = as.numeric(rma.ave[["beta"]]),
se = as.numeric(rma.ave[["se"]]),
ci.low = as.numeric(rma.ave[["ci.lb"]]),
ci.hi = as.numeric(rma.ave[["ci.ub"]])) )
# average for last 5 years only
dat.lim <- dat[dat$year>=2020,]
D.lim <- diag(dat.lim$SE)
R.lim <- R[12:16, 12:16]
V.lim <- D.lim %*% R.lim %*% D.lim
# intercept-only M.A: What is the pooled estimate across these 5 years?
rma.5yave <- rma.mv(yi = Beta,
V = V.lim,
data=dat.lim)
# save
group.B.5yave <- rbind(group.B.5yave, data.frame(
domain = domain,
group = group,
beta = as.numeric(rma.5yave[["beta"]]),
se = as.numeric(rma.5yave[["se"]]),
ci.low = as.numeric(rma.5yave[["ci.lb"]]),
ci.hi = as.numeric(rma.5yave[["ci.ub"]])) )
# include time as moderator to see if things have changed significantly for that group
# full time period
rma.trend.full <- rma.mv(yi = Beta,
V = V,
mods = ~ year.n,
data = dat)
# save
group.B.change <- rbind(group.B.change, data.frame(
domain=domain,
group=group,
beta=as.numeric(rma.trend.full[["beta"]][2]),
se=as.numeric(rma.trend.full[["se"]][2]),
ci.low=as.numeric(rma.trend.full[["ci.lb"]][2]),
ci.hi=as.numeric(rma.trend.full[["ci.ub"]][2]),
p=round(as.numeric(rma.trend.full[["pval"]][2]),4)))
}
}
# extract table of Ns per year
Ns <- data.frame(year = wgts.bygroup$year[wgts.bygroup$group=="Female"],
Females = wgts.bygroup$n[wgts.bygroup$group=="Female"],
Males = wgts.bygroup$n[wgts.bygroup$group=="Male"],
`10-12 years` = wgts.bygroup$n[wgts.bygroup$group=="10-12"],
`13-15 years` = wgts.bygroup$n[wgts.bygroup$group=="13-15"],
`High income` = wgts.bygroup$n[wgts.bygroup$group=="high.inc"],
`Middle income` = wgts.bygroup$n[wgts.bygroup$group=="mid.inc"],
`Low income` = wgts.bygroup$n[wgts.bygroup$group=="low.inc"]
)# approach: run one rma() per group factor (e.g. sex) and include an interaction between sex and time. Use custom function metareg.timeXgroup
# prepare the dataset of betas for the function
group.wgts.long$year.x <- group.wgts.long$year
group.wgts.long$group.factor <- factor(group.wgts.long$group, levels=c("10-12","13-15",
"Female" , "Male" ,
"high.inc", "mid.inc","low.inc"),
labels=c("age", "age", "sex", "sex", "income", "income", "income"))
group.wgts.long$dvlist <- factor(group.wgts.long$domain, levels=c("school.B","schoolwork.B","fam.B", "friend.B","appear.B"), labels=c("h.school.c", "h.schoolwork.c","h.fam.c", "h.friends.c", "h.appear.c"))
group.wgts.long$mean <- group.wgts.long$Beta
group.wgts.long$se <- group.wgts.long$SE
# limit to one group factor (CHECK metareg.timeXgroup() FUNCTION TO SEE WHY NOT WORKING (VAR-COVAR MATRIX))
df = group.wgts.long[group.wgts.long$group.factor=="sex",]
betas.sexdiffs <- metareg.timeXgroup(df, "full")
df = group.wgts.long[group.wgts.long$group.factor=="age",]
betas.agediffs <- metareg.timeXgroup(df, "full")
df = group.wgts.long[group.wgts.long$group.factor=="income",]
betas.incdiffs <- metareg.timeXgroup(df, "full")dplot <- group.wgts.long
dplot$group <- factor(dplot$group, levels=c("Female", "Male", "10-12", "13-15", "high.inc", "mid.inc","low.inc"))
ggplot(data=dplot[dplot$group=="Male" | dplot$group=="Female", ],
aes(x=year, y=Beta, group=group))+
geom_line(aes(colour=group), size=1)+
geom_ribbon(aes(ymin=Beta-2*SE, ymax=Beta+2*SE, fill=group),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
#scale_y_continuous(breaks=c(1,2), minor_breaks = c(0.5, 1.5,2.5),
# limits = c(0.3,2.6), name = "\nHappiness (-3 to +3)")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
#scale_shape_manual(values=c(3,15))+
scale_colour_manual(values=c("darkorange2", "turquoise3"))+
scale_fill_manual(values=c( "darkorange2", "turquoise3"))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=6))
ggplot(data=dplot[dplot$group=="10-12" | dplot$group=="13-15", ],
aes(x=year, y=Beta, group=group))+
geom_line(aes(colour=group), size=1)+
geom_ribbon(aes(ymin=Beta-2*SE, ymax=Beta+2*SE, fill=group),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
#scale_y_continuous(breaks=c(1,2), minor_breaks = c(0.5, 1.5,2.5),
# limits = c(0.3,2.6), name = "\nHappiness (-3 to +3)")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
#scale_shape_manual(values=c(3,15))+
scale_colour_manual(values=c("olivedrab3", "darkseagreen4"))+
scale_fill_manual(values=c("olivedrab3", "darkseagreen4"))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=6))
ggplot(data=dplot[str_detect(dplot$group, "inc"), ],
aes(x=year, y=Beta, group=group))+
geom_line(aes(colour=group), size=1)+
geom_ribbon(aes(ymin=Beta-2*SE, ymax=Beta+2*SE, fill=group),alpha=0.2)+
facet_grid(cols=vars(`Happiness with...`))+
theme_minimal()+
#scale_y_continuous(breaks=c(1,2), minor_breaks = c(0.5, 1.5,2.5),
# limits = c(0.3,2.6), name = "\nHappiness (-3 to +3)")+
scale_x_continuous(breaks=seq(2010,2024,2), limits = c(2009,2024), name = "")+
#scale_shape_manual(values=c(3,15))+
scale_colour_manual(values=c("darkslateblue","plum4","thistle" ))+
scale_fill_manual(values=c( "darkslateblue","plum4","thistle" ))+
theme(axis.text.x = element_text(angle=60, hjust = 1, size=6))
# averages over all 16 yrs
group.B.ave$domain.nice <- str_remove_all(group.B.ave$domain, ".B")
ggplot(data=group.B.ave[str_detect(group.B.ave$group, "-1"),],
aes(x= reorder(domain.nice,-beta), y=beta, ymin=ci.low, ymax=ci.hi, fill=domain.nice))+
geom_col()+
geom_errorbar(width=0)+
labs(x="Subdomains of life", y="Association with life satisfaction", title= "Split by age-group")+
facet_grid(~group)+
theme(axis.text.x = element_text(angle=30, size=9, hjust=1),
legend.position = "none")
ggplot(data=group.B.ave[str_detect(group.B.ave$group, "ale"),],
aes(x= reorder(domain.nice,-beta), y=beta, ymin=ci.low, ymax=ci.hi, fill=domain.nice))+
geom_col()+
geom_errorbar(width=0)+
labs(x="Subdomains of life", y="Association with life satisfaction", title= "Split by sex")+
facet_grid(~group)+
theme(axis.text.x = element_text(angle=30, size=9, hjust=1),
legend.position = "none")
ggplot(data=group.B.ave[str_detect(group.B.ave$group, "inc"),],
aes(x= reorder(domain.nice,-beta), y=beta, ymin=ci.low, ymax=ci.hi, fill=domain.nice))+
geom_col()+
geom_errorbar(width=0)+
labs(x="Subdomains of life", y="Association with life satisfaction", title= "Split by hsd income")+
facet_grid(~group)+
theme(axis.text.x = element_text(angle=30, size=9, hjust=1),
legend.position = "none")
# If we limit data to the last 5 years?
# average associational strengths in last 5 yrs
group.B.5yave$domain.nice <- str_remove_all(group.B.5yave$domain, ".B")
ggplot(data=group.B.5yave[str_detect(group.B.5yave$group, "-1"),],
aes(x= reorder(domain.nice,-beta), y=beta, ymin=ci.low, ymax=ci.hi, fill=domain.nice))+
geom_col()+
geom_errorbar(width=0)+
labs(x="Subdomains of life", y="Association with life satisfaction", title= "Split by age-group")+
facet_grid(~group)+
theme(axis.text.x = element_text(angle=30, size=9, hjust=1),
legend.position = "none")
ggplot(data=group.B.5yave[str_detect(group.B.5yave$group, "ale"),],
aes(x= reorder(domain.nice,-beta), y=beta, ymin=ci.low, ymax=ci.hi, fill=domain.nice))+
geom_col()+
geom_errorbar(width=0)+
labs(x="Subdomains of life", y="Association with life satisfaction", title= "Split by sex")+
facet_grid(~group)+
theme(axis.text.x = element_text(angle=30, size=9, hjust=1),
legend.position = "none")
ggplot(data=group.B.5yave[str_detect(group.B.5yave$group, "inc"),],
aes(x= reorder(domain.nice,-beta), y=beta, ymin=ci.low, ymax=ci.hi, fill=domain.nice))+
geom_col()+
geom_errorbar(width=0)+
labs(x="Subdomains of life", y="Association with life satisfaction", title= "Split by hsd income")+
facet_grid(~group)+
theme(axis.text.x = element_text(angle=30, size=9, hjust=1),
legend.position = "none")
There was no significant change over time in any group differences in associational strength between each sub-domain and general life satisfaction.
The association between subdomains and genral life happiness may be so high that there is not much room for moderation by group. But when we replace the life happiness item with SDQ total, we also don’t see much difference by sex.
Which groups and which domains showed a significant change over the 16 year time period?
ordered.df <- group.B.change[order(group.B.change$p),]
kable(ordered.df[1:10,], escape=F)%>%
kableExtra::kable_styling(full_width = F, font_size = 15)| domain | group | beta | se | ci.low | ci.hi | p | |
|---|---|---|---|---|---|---|---|
| 16 | appear.B | high.inc | 0.0099132 | 0.0026057 | 0.0048061 | 0.0150202 | 0.0001 |
| 6 | appear.B | 13-15 | 0.0085061 | 0.0024330 | 0.0037374 | 0.0132748 | 0.0005 |
| 11 | appear.B | Female | 0.0079901 | 0.0023101 | 0.0034624 | 0.0125177 | 0.0005 |
| 1 | appear.B | 10-12 | 0.0086877 | 0.0025798 | 0.0036314 | 0.0137439 | 0.0008 |
| 26 | appear.B | Male | 0.0079438 | 0.0027543 | 0.0025455 | 0.0133420 | 0.0039 |
| 21 | appear.B | low.inc | 0.0084731 | 0.0030817 | 0.0024330 | 0.0145131 | 0.0060 |
| 31 | appear.B | mid.inc | 0.0082887 | 0.0030383 | 0.0023337 | 0.0142436 | 0.0064 |
| 8 | friend.B | 13-15 | -0.0088178 | 0.0034441 | -0.0155681 | -0.0020675 | 0.0105 |
| 13 | friend.B | Female | -0.0084844 | 0.0033168 | -0.0149852 | -0.0019836 | 0.0105 |
| 12 | fam.B | Female | 0.0085772 | 0.0036735 | 0.0013772 | 0.0157772 | 0.0196 |
All groups showed a significant increase in the relative
“weighting” of appearance on general life satisfaction
(here’s the groups in order of that change):
- high income kids (B = 0.011)
- 13-15 year olds (B = 0.009)
- 10-12 year olds (B = 0.009)
- Females (B = 0.008)
- Males (B = 0.008)
- Low income kids (B = 0.009), p = .007
- Mid income kids (B = 0.007), p = .025
Partial & bivariate correlations (bounded between [-1, 1])
data4cor <- subset(usoc, select=c(h.life, h.appear, h.school, h.schoolwork, h.friends, h.fam))
data4cor <- data4cor[complete.cases(data4cor),]
names(data4cor) <- c("Life as whole", "Appearance", "School", "Schoolwork", "Friends", "Family")
# compute partial correlations
pc <- pcor(data4cor)
pc_matrix <- pc$estimate
# convert to edge list
edges <- pc_matrix %>%
as.data.frame() %>%
rownames_to_column("var1") %>%
pivot_longer(-var1, names_to = "var2", values_to = "partial_r") %>%
filter(var1 < var2)
# Define the other variables (in any order you prefer)
vars <- c("Life as whole", "Appearance", "School", "Schoolwork", "Friends", "Family")
others <- vars[vars != "Life as whole"]
# Create semicircle layout for the lower nodes
theta <- seq(-pi/6, -5*pi/6, length.out = length(others)) # wide, smooth arc
nodes <- data.frame(
name = vars,
x = c(0, # h.life centred in middle
cos(theta)), # other variables on arc
y = c(0.1, # h.life
-0.1 + 5*sin(theta)) # increase curvature for bottom row of nodes
)
g <- graph_from_data_frame(edges, directed = FALSE, vertices = nodes)
# Identify edges connected to "h.life"
connected_to_hlife <- incident(g, "Life as whole", mode = "all")
# Create an edge color vector for all edges
edge_cols <- rep("grey70", ecount(g))
edge_cols[connected_to_hlife] <- "black"
E(g)$edge_color <- edge_cols
ggraph(g, layout = "manual", x = x, y = y) +
geom_edge_link(aes(width = partial_r, colour= edge_color)) +
#geom_node_point(size = 12, colour = "white") +
geom_node_label(aes(label = name, fill=name), size =4, label.r = unit(7, "pt"),
label.padding = unit(6, "pt")) +
scale_edge_width(limits = c(0.01, 0.5), # partial correlation values
range = c(0.1, 5)) + # line thickness values
scale_fill_manual(values=c("#E78AC3FF" ,"#8DA0CBFF", "#FC8D62FF", "lightgrey","#66C2A5FF","#A6D854FF" ))+
scale_edge_colour_manual(values=c("black","grey80"))+
#scale_edge_colour_gradient2(low = "red", mid = "grey80", high = "blue") +
theme_void()+
theme(legend.position = "none",
plot.background = element_rect(fill="white"))
#print(plot)# run the network plot for every 5 years
plot_list <- vector("list", length = 3)
# for (YEAR in sort(unique(usoc$int.year.sched))){
for (window in sort(unique(usoc$time.window))) {
# select life satis data for 1 year only (no NAs)
data4cor <- subset(usoc[usoc$time.window==window,], select=c(h.life, h.appear, h.school, h.schoolwork, h.friends, h.fam))
data4cor <- data4cor[complete.cases(data4cor),]
# compute partial correlations
pc <- pcor(data4cor)
pc_matrix <- pc$estimate
# convert to edge list
edges <- pc_matrix %>%
as.data.frame() %>%
rownames_to_column("var1") %>%
pivot_longer(-var1, names_to = "var2", values_to = "partial_r") %>%
filter(var1 < var2)
# Define the other variables (in any order you prefer)
vars <- c("h.life", "h.appear" , "h.schoolwork","h.school", "h.friends", "h.fam" )
others <- vars[vars != "h.life"]
# Create semicircle layout for the lower nodes
theta <- seq(-pi/6, -5*pi/6, length.out = length(others)) # wide, smooth arc
nodes <- data.frame(
name = vars,
x = c(0, # h.life at centre top
cos(theta)), # other variables on arc
y = c(0.8, # h.life
-1 + 7*sin(theta)) # increase curvature for bottom row of nodes
)
g <- graph_from_data_frame(edges, directed = FALSE, vertices = nodes)
# Identify edges connected to "h.life"
connected_to_hlife <- incident(g, "h.life", mode = "all")
# Create an edge color vector for all edges
edge_cols <- rep("grey70", ecount(g))
edge_cols[connected_to_hlife] <- "black"
E(g)$edge_color <- edge_cols
plot_list[[as.numeric(window)]] <- ggraph(g, layout = "manual", x = x, y = y) +
geom_edge_link(aes(width = partial_r, colour= edge_color)) +
#geom_node_point(size = 12, colour = "white") +
geom_node_label(aes(label = name, fill=name), size =4, label.r = unit(7, "pt"),
label.padding = unit(6, "pt")) +
scale_edge_width(limits = c(0.05, 0.5), # partial correlation values
range = c(0.1, 4)) + # line thickness values
scale_fill_manual(values=c("#E78AC3FF" ,"#8DA0CBFF", "#FC8D62FF", "lightgrey","#66C2A5FF","#A6D854FF" ))+
scale_edge_colour_manual(values=c("black","grey80"))+
#scale_edge_colour_gradient2(low = "red", mid = "grey80", high = "blue") +
labs(title=paste("Time window:", window))+
theme_void()+
theme(legend.position = "none",
plot.background = element_rect(fill="white"))
#print(plot)
}
ggarrange(plot_list[[1]], plot_list[[2]], plot_list[[3]], nrow=1, ncol=3)
library(ppcor)
library(reshape2)
# save correlation between happiness in each domain & general life happiness
cortbl <- data.frame(year = sort(unique(usoc$int.year.sched)),
r.sch = NA,
r.schwk = NA,
r.fam = NA,
r.frnd = NA,
r.app = NA,
p.sch = NA,
p.schwk = NA,
p.fam = NA,
p.frnd = NA,
p.app = NA)
for (y in sort(unique(usoc$int.year.sched))){
# pearson's bivariate correlations
cortbl$r.sch[cortbl$year==y] <- cor(usoc$h.life[usoc$int.year.sched==y], usoc$h.school[usoc$int.year.sched==y],use="complete.obs")
cortbl$r.schwk[cortbl$year==y] <- cor(usoc$h.life[usoc$int.year.sched==y], usoc$h.schoolwork[usoc$int.year.sched==y],use="complete.obs")
cortbl$r.fam[cortbl$year==y] <- cor(usoc$h.life[usoc$int.year.sched==y], usoc$h.fam[usoc$int.year.sched==y],use="complete.obs")
cortbl$r.frnd[cortbl$year==y] <- cor(usoc$h.life[usoc$int.year.sched==y], usoc$h.friends[usoc$int.year.sched==y],use="complete.obs")
cortbl$r.app[cortbl$year==y] <- cor(usoc$h.life[usoc$int.year.sched==y], usoc$h.appear[usoc$int.year.sched==y],use="complete.obs")
# partial correlations
p <- pcor(na.omit(usoc[usoc$int.year.sched==y,
c("h.life", "h.school", "h.schoolwork", "h.fam", "h.friends", "h.appear")]))$estimate
p2 <- p
p2[lower.tri(p)] <- NA
p3 <- melt(p2, na.rm=T)
plife <- p3[p3$value!=1 & p3$Var1=="h.life",]
cortbl$p.sch[cortbl$year==y] <- plife$value[plife$Var2=="h.school"]
cortbl$p.schwk[cortbl$year==y] <- plife$value[plife$Var2=="h.schoolwork"]
cortbl$p.fam[cortbl$year==y] <- plife$value[plife$Var2=="h.fam"]
cortbl$p.frnd[cortbl$year==y] <- plife$value[plife$Var2=="h.friends"]
cortbl$p.app[cortbl$year==y] <- plife$value[plife$Var2=="h.appear"]
}
# plot these estimates as barchart
long.cor <- cortbl %>%
dplyr::select(year, r.sch, r.schwk ,r.fam ,r.frnd , r.app)%>%
gather(key = "domain", value="corr", -year)
long.partcor <- cortbl %>%
dplyr::select(year, p.sch, p.schwk ,p.fam ,p.frnd , p.app)%>%
gather(key = "domain", value="partialcorr", -year)
long.cor <- cbind(long.cor, long.partcor[,"partialcorr"])
names(long.cor)[str_detect(names(long.cor), "partialcorr")] <- "partialcorr"
long.cor$`Happiness with...` <- factor(long.cor$domain,
levels=c( "r.fam" ,"r.frnd" ,"r.sch", "r.schwk" , "r.app"),
labels=c( "Family" , "Friends",
"School", "Schoolwork",
"Appearance" ))
# bivariate correlations
ggplot(data=long.cor, aes(x=year, y=corr))+
geom_col(aes(fill=`Happiness with...`), alpha=0.7)+
geom_smooth(aes(x=year, y=corr), method="lm", formula=y~poly(x,2), colour="black", se=F)+
geom_smooth(aes(x=year, y=corr), colour="grey40", se=F)+
labs(y="Association with general life happiness (corr)\n", x="\nSurvey year", title="Bivariate correlations")+
facet_grid(~`Happiness with...`)+
#scale_y_continuous(limits=c(0, 0.5))+
scale_x_continuous(breaks= seq(2010,2024,2), limits = c(2008,2025), name = "")+
scale_fill_manual(values=cols5)+
scale_colour_manual(values=cols5)+
theme_minimal()+
theme(axis.text.x = element_text(size=7, angle=60, hjust=1, vjust = 1),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# partial correlations
ggplot(data=long.cor, aes(x=year, y=partialcorr))+
geom_col(aes(fill=`Happiness with...`), alpha=0.7)+
geom_smooth(aes(x=year, y=partialcorr), method="lm", formula=y~poly(x,2), colour="black", se=F)+
geom_smooth(aes(x=year, y=partialcorr), colour="grey40", se=F)+
labs(y="Association with general life happiness (pcorr)\n", x="\nSurvey year", title="Partial correlations")+
facet_grid(~`Happiness with...`)+
#scale_y_continuous(limits=c(0, 0.5))+
scale_x_continuous(breaks= seq(2010,2024,2), limits = c(2008,2025), name = "")+
scale_fill_manual(values=cols5)+
scale_colour_manual(values=cols5)+
theme_minimal()+
theme(axis.text.x = element_text(size=7, angle=60, hjust=1, vjust = 1),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# note: these partial correlations can't be compared directly with results of the svyglm
# plot network diagram for connection between domains
