Twine Manuscript Bayesian transition analyses
Mikael Rubin
January, 2021
A beginning section loads the necessary packages.
library(plyr)
library(dplyr)Load data
Data and models can be found on the Open Science Framework at https://osf.io/42ueg/.
transition1 <- read.csv("../DATA/transition_dat_t1.csv")
transition1$study <- "study1"
transition2 <- read.csv("../DATA/transition_dat_t2.csv")
transition2$study <- "study2"Create transition data frames
transitiondf <- plyr::rbind.fill(transition1, transition2)
transitiondf<-plyr::rename(transitiondf, c("Activities_1"="home_exercise",
"Activities_2"="out_exercise",
"Activities_3"="gym",
"Activities_4"="bake_cook",
"Activities_5"="volunteer",
"Activities_6"="television",
"Activities_7"="videogames",
"Activities_8"="knit",
"Activities_9"="read",
"Activities_10"="new_hobby",
"Activities_11"="alcohol",
"Activities_12"="marijuana",
"Activities_13"="other_substances",
"Activities_14"="reassurance",
"Activities_15"="social_media",
"Activities_16"="yoga",
"Activities_17"="writing",
"Activities_18"="play_music",
"Activities_19"="video_call_friends",
"Activities_20"="work_at_home",
"Activities_21"="work_outside_home"))
# Activities columns need to be 1,2,3,4 instead of 1,2,4,5
colnames(transitiondf)## [1] "X" "id" "EndDate"
## [4] "Age" "ethnicity" "race"
## [7] "ancestry" "sex" "sex_4_TEXT"
## [10] "gender" "work.impact" "social.impact"
## [13] "family.impact" "romantic.impact" "emotional.impact"
## [16] "health.impact" "home_exercise" "out_exercise"
## [19] "gym" "bake_cook" "volunteer"
## [22] "television" "videogames" "knit"
## [25] "read" "new_hobby" "alcohol"
## [28] "marijuana" "other_substances" "reassurance"
## [31] "social_media" "yoga" "writing"
## [34] "play_music" "video_call_friends" "work_at_home"
## [37] "work_outside_home" "PHQ.8_1" "PHQ.8_2"
## [40] "PHQ.8_3" "PHQ.8_4" "PHQ.8_5"
## [43] "PHQ.8_6" "PHQ.8_7" "PHQ.8_8"
## [46] "activities_activities" "consumable_activities" "tech_activities"
## [49] "activities_consumable" "consumable_consumable" "tech_consumable"
## [52] "activities_tech" "consumable_tech" "tech_tech"
## [55] "study"for(i in 27:47){
transitiondf[i] <- ifelse(transitiondf[i] == 4, c(3),
ifelse(transitiondf[i] == 5, c(4),
ifelse(transitiondf[i] == 1, c(1), c(2))))
}
#descriptive analyses:
demo <- transition1 %>% dplyr::select(Age:gender)
psych::describe(demo)## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf## vars n mean sd median trimmed mad min max range skew
## Age 1 262 18.98 1.41 19 18.71 1.48 18 30 12 3.07
## ethnicity 2 268 1.65 0.48 2 1.68 0.00 1 2 1 -0.61
## race* 3 268 11.18 4.50 13 11.62 2.97 1 15 14 -0.79
## ancestry* 4 266 3.92 7.36 1 1.75 0.00 1 31 30 2.38
## sex 5 268 1.29 0.46 1 1.24 0.00 1 2 1 0.91
## sex_4_TEXT 6 0 NaN NA NA NaN NA Inf -Inf -Inf NA
## gender* 7 267 12.06 8.99 9 11.27 10.38 1 34 33 0.70
## kurtosis se
## Age 16.08 0.09
## ethnicity -1.64 0.03
## race* -1.08 0.28
## ancestry* 4.13 0.45
## sex -1.17 0.03
## sex_4_TEXT NA NA
## gender* -0.48 0.55table(demo$race)##
## 1,4,6 1,6 2 2,3 2,3,6 2,5,6 2,6 3 3,4 3,4,6 3,6 4
## 1 1 1 61 1 1 1 6 16 1 2 3 54
## 4,6 6
## 18 101table(demo$ethnicity)##
## 1 2
## 95 173table(demo$gender)##
##
## 48
## bisexual
## 1
## Cis Male
## 1
## Cis-Woman
## 1
## cisgender
## 1
## Cisgender
## 4
## f
## 1
## female
## 61
## Female
## 52
## Female
## 6
## Female, although I don't really care if you use they/them or call me bro or dude.
## 1
## feminine
## 1
## He/Him
## 1
## He/Him. Male.
## 1
## Her
## 1
## her/she
## 1
## Hetero Male
## 1
## heterosexual
## 2
## Heterosexual
## 3
## Heterosexual Male
## 1
## male
## 16
## Male
## 37
## Male
## 2
## Male (Heterosexual)
## 1
## man
## 1
## Man
## 2
## Non binary
## 1
## Non-binary
## 1
## she/her/hers
## 1
## Straight
## 1
## Straight male
## 1
## woman
## 5
## Woman
## 8
## Woman
## 1table(demo$sex)##
## 1 2
## 190 78demo1 <- transition1 %>% dplyr::select(Age:gender)
table(demo1$race)##
## 1,4,6 1,6 2 2,3 2,3,6 2,5,6 2,6 3 3,4 3,4,6 3,6 4
## 1 1 1 61 1 1 1 6 16 1 2 3 54
## 4,6 6
## 18 101table(demo1$ethnicity)##
## 1 2
## 95 173table(demo1$gender)##
##
## 48
## bisexual
## 1
## Cis Male
## 1
## Cis-Woman
## 1
## cisgender
## 1
## Cisgender
## 4
## f
## 1
## female
## 61
## Female
## 52
## Female
## 6
## Female, although I don't really care if you use they/them or call me bro or dude.
## 1
## feminine
## 1
## He/Him
## 1
## He/Him. Male.
## 1
## Her
## 1
## her/she
## 1
## Hetero Male
## 1
## heterosexual
## 2
## Heterosexual
## 3
## Heterosexual Male
## 1
## male
## 16
## Male
## 37
## Male
## 2
## Male (Heterosexual)
## 1
## man
## 1
## Man
## 2
## Non binary
## 1
## Non-binary
## 1
## she/her/hers
## 1
## Straight
## 1
## Straight male
## 1
## woman
## 5
## Woman
## 8
## Woman
## 1table(demo1$sex)##
## 1 2
## 190 78game.feedback <- transition1 %>% dplyr::select(PHQ.8_1:PHQ.8_8)
psych::alpha(game.feedback)##
## Reliability analysis
## Call: psych::alpha(x = game.feedback)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.89 0.89 0.89 0.5 8 0.01 0.87 0.68 0.5
##
## lower alpha upper 95% confidence boundaries
## 0.87 0.89 0.91
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## PHQ.8_1 0.87 0.87 0.86 0.49 6.7 0.012 0.0069 0.49
## PHQ.8_2 0.86 0.87 0.86 0.48 6.5 0.012 0.0058 0.48
## PHQ.8_3 0.88 0.88 0.87 0.52 7.5 0.011 0.0069 0.51
## PHQ.8_4 0.87 0.87 0.86 0.49 6.6 0.012 0.0083 0.49
## PHQ.8_5 0.87 0.88 0.87 0.51 7.2 0.011 0.0091 0.51
## PHQ.8_6 0.87 0.87 0.86 0.49 6.7 0.012 0.0082 0.49
## PHQ.8_7 0.88 0.88 0.87 0.51 7.3 0.011 0.0084 0.49
## PHQ.8_8 0.88 0.89 0.88 0.53 7.8 0.011 0.0068 0.53
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PHQ.8_1 268 0.79 0.79 0.76 0.71 0.81 0.84
## PHQ.8_2 268 0.82 0.82 0.81 0.76 0.85 0.87
## PHQ.8_3 268 0.71 0.69 0.63 0.60 1.05 1.03
## PHQ.8_4 268 0.80 0.80 0.77 0.73 1.21 0.92
## PHQ.8_5 268 0.74 0.73 0.68 0.64 1.04 1.02
## PHQ.8_6 268 0.80 0.79 0.76 0.72 0.90 0.97
## PHQ.8_7 268 0.71 0.72 0.66 0.61 0.82 0.96
## PHQ.8_8 268 0.62 0.66 0.58 0.55 0.29 0.62
##
## Non missing response frequency for each item
## 0 1 2 3 miss
## PHQ.8_1 0.41 0.40 0.14 0.04 0
## PHQ.8_2 0.41 0.40 0.13 0.06 0
## PHQ.8_3 0.37 0.33 0.16 0.13 0
## PHQ.8_4 0.22 0.46 0.19 0.12 0
## PHQ.8_5 0.37 0.34 0.16 0.13 0
## PHQ.8_6 0.43 0.32 0.15 0.09 0
## PHQ.8_7 0.47 0.32 0.12 0.09 0
## PHQ.8_8 0.77 0.19 0.02 0.02 0psych::describe(game.feedback)## vars n mean sd median trimmed mad min max range skew kurtosis se
## PHQ.8_1 1 268 0.81 0.84 1 0.71 1.48 0 3 3 0.82 0.03 0.05
## PHQ.8_2 2 268 0.85 0.87 1 0.74 1.48 0 3 3 0.84 0.00 0.05
## PHQ.8_3 3 268 1.05 1.03 1 0.94 1.48 0 3 3 0.61 -0.81 0.06
## PHQ.8_4 4 268 1.21 0.92 1 1.14 1.48 0 3 3 0.48 -0.56 0.06
## PHQ.8_5 5 268 1.04 1.02 1 0.93 1.48 0 3 3 0.64 -0.73 0.06
## PHQ.8_6 6 268 0.90 0.97 1 0.76 1.48 0 3 3 0.80 -0.44 0.06
## PHQ.8_7 7 268 0.82 0.96 1 0.67 1.48 0 3 3 0.97 -0.09 0.06
## PHQ.8_8 8 268 0.29 0.62 0 0.16 0.00 0 3 3 2.51 6.85 0.04game.feedback2 <- transition2 %>% dplyr::select(PHQ.8_1:PHQ.8_8)
psych::alpha(game.feedback2)##
## Reliability analysis
## Call: psych::alpha(x = game.feedback2)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.89 0.89 0.89 0.5 8.1 0.011 0.97 0.7 0.51
##
## lower alpha upper 95% confidence boundaries
## 0.87 0.89 0.91
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## PHQ.8_1 0.88 0.88 0.87 0.51 7.4 0.012 0.0049 0.51
## PHQ.8_2 0.87 0.87 0.86 0.49 6.7 0.013 0.0046 0.50
## PHQ.8_3 0.88 0.88 0.87 0.50 7.1 0.012 0.0053 0.51
## PHQ.8_4 0.87 0.87 0.86 0.49 6.8 0.013 0.0046 0.50
## PHQ.8_5 0.88 0.88 0.87 0.50 7.1 0.012 0.0072 0.51
## PHQ.8_6 0.87 0.88 0.87 0.50 7.0 0.012 0.0058 0.50
## PHQ.8_7 0.87 0.87 0.87 0.50 7.0 0.012 0.0070 0.50
## PHQ.8_8 0.89 0.89 0.88 0.53 7.9 0.011 0.0032 0.51
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PHQ.8_1 229 0.71 0.72 0.67 0.63 0.97 0.83
## PHQ.8_2 229 0.80 0.80 0.78 0.73 0.91 0.85
## PHQ.8_3 229 0.76 0.75 0.71 0.67 1.13 0.98
## PHQ.8_4 229 0.80 0.79 0.77 0.72 1.40 0.98
## PHQ.8_5 229 0.76 0.75 0.70 0.67 1.02 1.00
## PHQ.8_6 229 0.77 0.77 0.73 0.68 0.97 0.99
## PHQ.8_7 229 0.78 0.77 0.73 0.69 1.01 1.01
## PHQ.8_8 229 0.64 0.66 0.58 0.55 0.36 0.72
##
## Non missing response frequency for each item
## 0 1 2 3 miss
## PHQ.8_1 0.30 0.49 0.15 0.06 0
## PHQ.8_2 0.36 0.42 0.17 0.05 0
## PHQ.8_3 0.31 0.36 0.21 0.11 0
## PHQ.8_4 0.18 0.41 0.23 0.17 0
## PHQ.8_5 0.38 0.35 0.16 0.12 0
## PHQ.8_6 0.41 0.32 0.18 0.10 0
## PHQ.8_7 0.38 0.34 0.16 0.12 0
## PHQ.8_8 0.75 0.18 0.04 0.03 0psych::describe(game.feedback2)## vars n mean sd median trimmed mad min max range skew kurtosis se
## PHQ.8_1 1 229 0.97 0.83 1 0.89 1.48 0 3 3 0.69 0.04 0.06
## PHQ.8_2 2 229 0.91 0.85 1 0.83 1.48 0 3 3 0.64 -0.29 0.06
## PHQ.8_3 3 229 1.13 0.98 1 1.04 1.48 0 3 3 0.46 -0.84 0.06
## PHQ.8_4 4 229 1.40 0.98 1 1.37 1.48 0 3 3 0.26 -0.95 0.06
## PHQ.8_5 5 229 1.02 1.00 1 0.90 1.48 0 3 3 0.66 -0.67 0.07
## PHQ.8_6 6 229 0.97 0.99 1 0.84 1.48 0 3 3 0.67 -0.67 0.07
## PHQ.8_7 7 229 1.01 1.01 1 0.89 1.48 0 3 3 0.67 -0.68 0.07
## PHQ.8_8 8 229 0.36 0.72 0 0.19 0.00 0 3 3 2.19 4.41 0.05Bayesian regression
transitiondf$study.f <- as.factor(transitiondf$study)
transitiondf$id.f <- as.factor(transitiondf$id)
Prior <- brms::set_prior("normal(0,1)", class= "b")Fit Bayesian GLMMMs
Note that this will take some time—as such, we include the saved models. If running this code directly, you may run each code chunk below, or set create_models to true in the preamble at the top of this document.
load("../DATA/models/brm_models.Rdata")out1 <- brms::brm(tech_tech ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out2 <- brms::brm(activities_activities ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out3 <- brms::brm(consumable_activities ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out4 <- brms::brm(tech_activities ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out5 <- brms::brm(activities_consumable ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out6 <- brms::brm(consumable_consumable ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out7 <- brms::brm(tech_consumable ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out8 <- brms::brm(activities_tech ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out9 <- brms::brm(consumable_tech ~ 0 + study.f, data=transitiondf, prior = Prior, sample_prior = TRUE)
out10 <- brms::brm(phq_tot ~ tech_tech + activities_activities + consumable_activities+
tech_activities + activities_consumable + consumable_consumable + tech_consumable +
activities_tech + consumable_tech + study.f,
data=transitiondf, seed = 1, iter=8000){BRMS} hypothesis testing
h1 <- brms::hypothesis(out1, "study.fstudy1 = study.fstudy2")
h2 <- brms::hypothesis(out2, "study.fstudy1 = study.fstudy2")
h3 <- brms::hypothesis(out3, "study.fstudy1 = study.fstudy2")
h4 <- brms::hypothesis(out4, "study.fstudy1 = study.fstudy2")
h5 <- brms::hypothesis(out5, "study.fstudy1 = study.fstudy2")
h6 <- brms::hypothesis(out6, "study.fstudy1 = study.fstudy2")
h7 <- brms::hypothesis(out7, "study.fstudy1 = study.fstudy2")
h8 <- brms::hypothesis(out8, "study.fstudy1 = study.fstudy2")
h9 <- brms::hypothesis(out9, "study.fstudy1 = study.fstudy2")
1/h1$hypothesis$Evid.Ratio## [1] 0.013242511/h2$hypothesis$Evid.Ratio## [1] 0.020904471/h3$hypothesis$Evid.Ratio## [1] 0.0058455951/h4$hypothesis$Evid.Ratio## [1] 0.003113781/h5$hypothesis$Evid.Ratio## [1] 0.0030620711/h6$hypothesis$Evid.Ratio## [1] 0.0035836951/h7$hypothesis$Evid.Ratio## [1] 0.002924771/h8$hypothesis$Evid.Ratio## [1] 0.0028936081/h9$hypothesis$Evid.Ratio## [1] 0.03810143Plot
out1.loo<-brms::loo(out1)
brms::pp_check(out1)## Using 10 posterior samples for ppc type 'dens_overlay' by default.
out1.loo<-brms::loo(out10)
brms::pp_check(out10)## Using 10 posterior samples for ppc type 'dens_overlay' by default.