Main analysis
1. Time trends of life satisfaction
Extract means & proportions from each calendar year (intercept-only models)
# just for satisfaction with life as a whole, create a plot per year and ADD THE MEAN AS A VERTICAL REF LINE
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"))
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/lifesatis_distribution_change.tiff", units="in", width=8, height=6, dpi=600, compression = 'lzw')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.
The GLS-AR(1) model says: * “Once I accept that this population evolves smoothly over time, there is much less independent information than it looks like.” * “Under a model that assumes strong, persistent dependence across all calendar years, the data do not contain enough independent information to estimate the slope precisely.” * That model assumption is not aligned with your sampling design, so its inferential pessimism is not necessarily meaningful.
Notes on time-dependency
At the individual level: - Each person contributes up to 6 consecutive years. - Entry and exit are staggered across calendar time.
At the year-level aggregate (your observed means): - Adjacent years share many individuals. - Years 5–6 apart still share some individuals. - Very distant years share few or no individuals. - The overlap decays approximately linearly, not geometrically.
Crucially: - Dependence is finite-range and smooth - There is no single correlation parameter that governs all lags - Long-lag correlations are weak but non-zero in principle
This is not an AR(1) process.
# TEST SIGNIFICANCE OF LINEAR TRENDS
# 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 M.A: 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...
# 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## Warning: 1 row with NAs omitted from model fitting.
## Warning: 1 row with NAs omitted from model fitting.
## Warning: 1 row with NAs omitted from model fitting.
## Warning: 1 row with NAs omitted from model fitting.
## Warning: 1 row with NAs omitted from model fitting.
## Warning: 1 row with NAs omitted from model fitting.# 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)))
}mean(prop.fullsample$prop[prop.fullsample$domain=="unhappy.life" & (prop.fullsample$windows==1)], na.rm=T)## [1] 0.03932591mean(prop.fullsample$prop[prop.fullsample$domain=="unhappy.life" & (prop.fullsample$windows==2)], na.rm=T)## [1] 0.05639085mean(prop.fullsample$prop[prop.fullsample$domain=="unhappy.life" & (prop.fullsample$windows==3)], na.rm=T)## [1] 0.07472576# TEST SIGNIFICANCE OF LINEAR TRENDS
# 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
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))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/fullpop_trends.tiff", units="in", width=6, height=5, dpi=600, compression = 'lzw')
# 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))## `geom_smooth()` using formula = 'y ~ x'
# ...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()
# MEDIANS
# means.fullsample$median[means.fullsample$domain=="h.appear.M"] = means.fullsample$median[means.fullsample$domain=="h.appear.M"]+0.05
#
# means.fullsample$median[means.fullsample$domain=="h.schoolwork.M"] = means.fullsample$median[means.fullsample$domain=="h.schoolwork.M"]-0.05
#
# means.fullsample$median[means.fullsample$domain=="h.fam.M"] = means.fullsample$median[means.fullsample$domain=="h.fam.M"]+0.05
#
# ggplot(data=means.fullsample, aes(x=year.x, y=median, group=domain, fill=domain))+
# geom_point(aes(colour=domain))+
# geom_line(aes(colour=domain), size=1.5, alpha=0.8)+
# labs(y="Median happiness (-3 to +3)")+
# lims(y=c(0, 3.5))+
# theme(axis.text.x = element_text(angle=30, hjust = 1))Note wgtd median vs means - The weighted mean averages across the scale and can take on decimals (e.g., 1.42), even though no single respondent chose that exact value. - The weighted median, by contrast, has to land on a real category, because it’s defined as the “middle observation” once all weights are accounted for.
Precision It seems my mean estimates are becoming less precise over time. This is 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.
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))
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/fullpop_Variance.tiff", units="in", width=8, height=5, dpi=600, compression = 'lzw')Note that the y-axis units here have no meaningful value, they are just in terms of the likert satisfaction scale
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))
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/fullpop_unhappytrends_line.tiff", units="in", width=8, height=5, dpi=600, compression = 'lzw')
# 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")
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/fullpop_unhappytrends_bars_rows.tiff", units="in", width=3, height=8, dpi=600, compression = 'lzw')
# 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")
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/fullpop_unhappytrends_bars_grid.tiff", units="in", width=6, height=5, dpi=600, compression = 'lzw')# 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")
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/trends_bygroup.tiff", units="in", width=9, height=8, dpi=600, compression = 'lzw')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))
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/alltrends_by3groupage.tiff", units="in", width=12, height=5, dpi=600, compression = 'lzw')
# "tomato", "tomato4","cyan2", "olivedrab2", "darkseagreen3","aquamarine4"
# ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/alltrends_byspecif_age.tiff", units="in", width=12, height=5, dpi=600, compression = 'lzw')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))
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/trends_bycountry.tiff", units="in", width=12, height=5, dpi=600, compression = 'lzw')
# sample sizes are too small-- trends are v noisy
# calculate 2-year averages / extract means per wave.
# and re-run plots
#comparing scotland and england on school vars might be a bit more reliable
# what's NI doing to make their schools more likable?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())
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/sexage_trends_4plot.tiff", units="in", width=11, height=5, dpi=600, compression = 'lzw')
# 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
Beta coefficients (-Inf, Inf)
– full population
# 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 (univariate-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
# points connected by straight lines
# ggplot(data=wgts.long, aes(x=year, y=Beta, group=`Happiness with...`, fill=`Happiness with...`))+
# geom_ribbon(aes(ymin=Beta-SE, ymax=Beta+SE), alpha=0.5)+
# geom_point(aes(colour=`Happiness with...`))+
# geom_line(aes(colour=`Happiness with...`))+
# labs(y="Association with gen life happiness (B ± SE)")+
# scale_colour_manual(values=cols5)+
# scale_fill_manual(values=cols5)+
# theme_minimal()+
# theme(axis.text.x = element_text(angle=30, hjust = 1))
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())
ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/weights_over_time_unadj.tiff", units="in", width=10, height=5, dpi=600, compression = 'lzw')
# 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())
#ggsave("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/weights_over_time.tiff", units="in", width=10, height=5, dpi=600, compression = 'lzw')
## 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
But these correlations are relative on each other… what is the pure bivariate correlation between each domain and general life satisfaction
# 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)library(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)))
}
}
save(group.wgts.long, group.B.ave, group.B.5yave, group.B.change,
file="/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/relat.wgts.groups.Rdata")
# 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"]
)load("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/relat.wgts.groups.Rdata")# a different approach: run one rma() per group factor (e.g. sex) and include an interaction between sex and time
# see if your function metareg.timeXgroup() can be used for this...
# 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")
# pretty sure reference groups are: females, 10-12 yo, high income -- see plots for proofAdapt summary below The plots below show that there are no strong group differences in B-estimates (how sub-domains relate to general life happiness), nor in the change in those association over time (for example, appearance B increases in all groups)
The relative importance of family was the only domain where there was a significant (p<.05) moderation. Family had weaker associations with life satisfaction for males compared to females (but p = .04) and for older compared to younger kids (p = .01). We do not see the same patterns when using SDQ tot as the outcome.
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))
# make middle income go in the middle!!
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))
Further take-aways: - these association estimates are very
unstable. They fluctuate greatly from year to year, which makes any
group difference unlikely to be detected. - but all we can conclude from
this data is that there are no significant differences (that can be
detected) in the influence of each area of life on general life
happiness. - the beta estimates for each subdomain are so remarkably
similar across groups, that it actyally makes me think this is a
statistical similarity rather than actual “influence” of each domain.
i.e. school & schoolwork will be simmilarly ranked, so they are
downweighted in tandem. Family & appearance are the most
“independent” or different constructs of the 5, so they “show up”. Try
and phrase this in a more scientific way for the paper.
But are we more interested in the order of aspects of life, for different subgroups?
# 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")
When we had statistically looked at the beta estimate for group, there was no significant change over time in any group differences in associational strength between each sub-domain and general life satisfaction. Sensitivity analysis with SDQ only showed 1 change over time: the association between family and life happiness growing in the older kids more than the younger kids (but only p = .03)
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
Then there were borderline sig decreases in relative association between friends & general life happiness: - females (B = -0.007) - 13-15 year olds (B = -0.009)
Partial & bivariate correlations (-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)
ggsave(paste0("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/networkplots/avg_fullperiod.tiff"), units="in", width=8, height=7, dpi=600, compression = 'lzw')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)
ggsave(paste0("/Users/niamhdooley/Library/CloudStorage/OneDrive-RoyalCollegeofSurgeonsinIreland/DOROTHY CO-FUND/WP2-Understanding Society/Project A. Targets for intervention in life satisfaction/R outputs/networkplots/3windows.tiff"), units="in", width=20, height=5, dpi=600, compression = 'lzw')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 domainsYou need to think about what this means: - the bivariate correlation between all subdomains and happiness with life in general is increasing over time…. these 5 domain have become more important? What hasnt? Or all are subject to higher scores (i.e. satisfacation with everything is decreasing) -
We can see that the % of minority ethnicities “peak” around 2015 and
2023 around the times of IEM sample boosting. The weights seem to adjust
for that design-related bias. Even with the weighting, the proportion of
ethnic minorities (all others not British/Irish) rise over time, and
rise steeply from 2021.
