Data and models for NRSA power analysis

Study description

Data for this project come from a larger parent study (“The effect of intranasal oxytocin on neural functioning in widow(er)s” (#1503744420, PI: Mary-Frances O’Connor, Co-PI: Brian Arizmendi). It was a double-blind placebo-controlled randomized crossover study, designed to test the neuropeptide oxytocin as a potential mechanism of poor adaptation in complicated vs. typical grief. Participants were 40 bereaved adults ages 55-80, who experienced the death of a spouse or long-term romantic partner between six months and three years prior to their participation.

Participants attended two identical experimental sessions, at which they received one of two intranasal sprays (treatment A/B) and then took part in an fMRI scan, which involved a structural sequence, a behavioral task involving viewing images of their deceased partner, and a resting state sequence (in that order).

Data wrangling steps

Import static and dynamic FNC from GIFT output

library(tidyverse)
filter <- dplyr::filter
select <- dplyr::select

# read in static FNC stats
## FNC = correlation between pair of IC
fnc_corrs_all <- read_csv("~/Dropbox/fnc_corrs_all.csv", col_names = TRUE)

# subset only FNCs of interest by session

# placebo session only:
fnc_corrs_all_Pl <- select(fnc_corrs_all,c(ID, 
fnc_c04c13_txA, # intra-frontoparietal network (R, L dlPFC) / Placebo
fnc_c04c14_txA, # R frontoparietal-default mode network / Placebo
fnc_c14c17_txA, # intra-default mode network / Placebo
fnc_c14c27_txA, # intra-default mode network / Placebo
fnc_c19c27_txA)) # salience network-default mode network / Placebo

# drop the "txA" suffix and add a new column "tx" for the long dataframe
fnc_corrs_all_Pl <- fnc_corrs_all_Pl %>%
   rename_all(.funs = funs(sub("\\_txA", "", names(fnc_corrs_all_Pl))))
fnc_corrs_all_Pl$tx <- "txA"
names(fnc_corrs_all_Pl)

# oxytocin session only:
fnc_corrs_all_Ot <- select(fnc_corrs_all,c(ID, 
fnc_c04c13_txB, # intra-frontoparietal network (R, L dlPFC) / Oxytocin
fnc_c04c14_txB, # R frontoparietal-default mode network / Oxytocin
fnc_c14c17_txB, # intra-default mode network / Oxytocin
fnc_c14c27_txB, # intra-default mode network / Oxytocin
fnc_c19c27_txB)) # salience network-default mode network / Oxytocin

# drop the "txB" suffix and add a new column "tx" for the long dataframe
fnc_corrs_all_Ot <- fnc_corrs_all_Ot %>%
   rename_all(.funs = funs(sub("\\_txB", "", names(fnc_corrs_all_Ot))))
fnc_corrs_all_Ot$tx <- "txB"
names(fnc_corrs_all_Ot)

# create the long FNC dataframe
fnc_corrs_all_long <- bind_rows(fnc_corrs_all_Pl, fnc_corrs_all_Ot)

# read in the temporal dynamic functional connectivity (DFNC) stats
# already long so no reshaping needed
state_vector_stats <- read_csv("~/Dropbox/state_vector_stats.csv", col_names = TRUE)
head(state_vector_stats)

# merge on columns ID and tx
gift_out <- left_join(state_vector_stats,fnc_corrs_all_long, by=c("ID", "tx"))

Add state and time-invariant variables

# using the % of content that was spouse/death-related from their post scan reports as a very imperfect measure of ruminative content for Aim 1 H3
# (this will probably change)

# placebo:
rumA <- read_csv("/Users/sarenseeley/Dropbox/Grants\ and\ Applications/NRSA\ F31\ 2018-2019/Resubmission\ Reviews/power-analysis/rumination_txA.csv")

# create composite variables for (1) focus on spouse in post-scan responses, and (2) of all of the things they reported thinking about, what % of what they mentioned was related to their spouse/death?
rumA$focus_spouse <- (rumA$spouse_1+rumA$spouse_2) 
rumA$focus_task <- (rumA$task_1+rumA$task_2) 
rumA$focus_whoto <- (rumA$whoto_1+rumA$whoto_2) 
rumA$focus_otherpics <- (rumA$otherpics_1+rumA$otherpics_2) 
rumA$focus_nonspouse <- (rumA$spouse_1+rumA$spouse_2+rumA$task_1+rumA$task_2+rumA$whoto_1+rumA$whoto_2+rumA$otherpics_1+rumA$otherpics_2+rumA$somatic_1+rumA$somatic_2+rumA$metacog_1+rumA$metacog_2)
rumA$focus_spouse_perc <- (rumA$focus_spouse/rumA$focus_nonspouse)
View(rumA)

# replace NaNs (generated by dividing zero by zero) with 0
fix_nan <- function(x){
    x[is.nan(x)] <- 0
    x
}
rumA$focus_spouse_perc <- fix_nan(rumA$focus_spouse_perc)   

# drop variables extraneous for power analysis
rumA <- rumA %>% select(-c(post_sroe_1_A, spouse_1, task_1, whoto_1, otherpics_1, somatic_1, metacog_1, post_sroe_2_A, spouse_2, task_2, whoto_2, otherpics_2, somatic_2, metacog_2, focus_spouse, focus_task, focus_whoto, focus_otherpics, focus_nonspouse))
names(rumA)

                                
# make a subset of data from master dataset
dataforjf <- select(data_spr19,c(ID, 
                                 group, # complicated grief (CG) or non-complicated grief (NCG), based on tot_icg clinical cutoff of >25
                                 tot_icg, # Inventory of Complicated Grief total score
                                 tot_bdi, # Beck Depression Inventory total
                                 age_yrs, # age in years
                                 sex_m, # sex (m is reference category)
                                 timesincedeath, # time elapsed (in days) between spouse's death and study participation
                                 tx_v1)) # which treatment they got at their first visit (A = placebo)

# replace "D" prefix in ID to "sub-" so that the ID variables match across dataframes to merge
dataforjf <- dataforjf %>% 
  mutate(ID = str_replace(ID, "D", "sub-"))
head(dataforjf)

# merge state reports with FNC/DFNC variables
scandata <- left_join(gift_out, rumA, by=c("ID","tx"))

# merge with time-invariant variables
dataforjf <- left_join(dataforjf,scandata, by=c("ID"))

# recode factors: substitute "placebo" for A/txA, and "oxytocin" for B/txB
# placebo = 0, oxytocin = 1
dataforjf <- dataforjf %>% mutate(tx_v1 = recode_factor(tx_v1, "A" = 0, "B" = 1), tx = recode_factor(tx, "txA" = 0, "txB" = 1)) 

# save as csv file
write_csv(dataforjf, "/Users/sarenseeley/Dropbox/Grants\ and\ Applications/NRSA\ F31\ 2018-2019/Resubmission\ Reviews/power-analysis/dataforjf.csv")

Which states do I care about?

Looking at results from 1-sample t test (FDR-corrected p < .05) in tDFNC toolbox…

State 1 is characterized by:

• greater + connectivity between default mode (IC14, IC17 [precuneus/PCC, autobio memory])
• greater + connectivity between the salience network (IC19) and right frontoparietal network (IC4)
• greater + connectivity between the salience network (IC19) and default mode network (IC14 [precuneus/PCC])
• greater + connectivity between the salience network (IC19) and default mode network (IC14 [precuneus/PCC])
• greater - connectivity between the salience network (IC19) and default mode network (IC27 [midline/self-reference/social cog])

State 2 is characterized by strong intercorrelations between and within all networks except for:

• between the salience network (IC19) and one of the default mode ICs (IC27 [midline/self-reference/social cog])
• IC27 is also less strongly connected with the right frontoparietal network (IC13) and another default mode IC (IC17 [autobio memory])

 

State 3 is characterized by:

• greater + connectivity within the default mode network (ICs 14, 17, 27)
• greater - connectivity between the salience network (IC19) and default mode network (IC27 [midline/self-reference/social cog])

State 4 is characterized by:

• greater + connectivity within the frontoparietal network (ICs 4, 13)
• greater + connectivity within the default mode ICs (IC14 [precuneus/PCC] w/IC27 [midline/self-reference/social cog] and IC14 [precuneus/PCC] w/IC17 [autobio memory])
• greater + connectivity between the frontoparietal network (IC 4, 13) and midline default mode network (IC14 [precuneus/PCC], IC27 [midline/self-reference/social cog])
• greater - connectivity between the salience network (IC19) and midline default mode network (IC14 [precuneus/PCC], IC27 [midline/self-reference/social cog])
• greater - connectivity between the salience network (IC19) and left frontoparietal network (IC4)

Interpretations

State 1: Vast majority (35/38) of subjects show this state. This one seems mostly about connectivity within the posterior DMN, and of L frontoparietal with default mode and salience networks.

State 2: Strong intercorrelations within and between nearly all ICs. (Meaning???) This is the state in which default mode components (IC14, IC17) are positively coupled with the salience network component (IC19), and looks the most like Fig. 2 from my proposal:

Figure 2. Hypothesized major network interactions in complicated grief, showing the flow and content of spontaneous thought as potentially more automatically constrained and less variable due to (a) lesser influence from executive control regions over default mode influence26, (b) greater influence from default mode core subsystem regions48, and (c) greater influence of the salience network – so thoughts are intrusive, inflexible, and fixated on the loss.

State 3: Cognition driven by internally-oriented thoughts from default mode network (memory, self-referential). Salience and frontoparietal networks are negatively correlated with default mode activity - few constraints over thought. Mind-wandering?

State 4: Strong FNC within and between default mode and frontoparietal control network components. Salience network negatively coupled with frontoparietal network and default mode components. May reflect more goal-directed cognition?

Variables in dataset ‘dataforjf.csv’

ID: Subject ID

tot_icg: Inventory of Complicated Grief total score (grief symptoms measure)

group: Complicated grief (CG) or non-complicated grief (NCG), based on tot_icg clinical cutoff of >25 (NCG = 0, CG = 1)

age_yrs: Age in years

timesincedeath: Time (in days) elapsed between date of spouse’s death and date of study participation (session 1)

sex_m: Sex (male = 0, female = 1)

tx_v1: Treatment received at first visit (placebo = 0, oxytocin = 1) (for order effects)

tx: Observation comes from which session? (placebo = 0, oxytocin = 1)

focus_spouse_perc: Of the number of different content areas mentioned in the post-scan free response self-report questions about what they were thinking/feeling and doing during the task, what proportion was spouse-related? (each of six content areas was scored 0/1 for absent/present, scores on spouse dummy variable summed [range 0-2], then divided by an aggregate score that summed all dummy variables to yield a measure of how focused they were on their spouse and/or their death during the task.) I need to think more about how best to get a metric that captures what I’m looking for – I’m not sure this is it.

meta_states_numstates: Number of unique windows for each subject

meta_states_totaldist: Sum of L1 distances between successive meta-states for each subject

meta_states_changestates: Number of times each subject changes from one meta state to other

meta_states_statespan: Maximum L1 distance between states for each subject

statestats_fract_time_s1: Fraction of time spent in state 1

statestats_fract_time_s2: Fraction of time spent in state 2

statestats_fract_time_s3: Fraction of time spent in state 3

statestats_fract_time_s4: Fraction of time spent in state 4

statestats_dwelltime_s1: The average amount of time spent state 1 before shifting to another state (for each subject and session)

statestats_dwelltime_s2: The average amount of time spent state 2 before shifting to another state (for each subject and session)

statestats_dwelltime_s3: The average amount of time spent state 3 before shifting to another state (for each subject and session)

statestats_dwelltime_s4: The average amount of time spent state 4 before shifting to another state (for each subject and session)

statestats_numtransitions: How many state transitions occurred (for each subject and session)

fnc_c04c13: Static FNC within the frontoparietal network

fnc_c04c14: Static FNC between the frontoparietal and default mode network

fnc_c14c17: Static FNC within the default mode network

fnc_c14c27: Static FNC within the default mode network

fnc_c19c27: Static FNC between the salience and default mode network

Aims, analyses and syntax

Aim 1

Aim 1. To test whether dynamic functional connectivity under placebo differs among widowed older adults who are adjusting well, versus those who are adjusting poorly (i.e., complicated grief).

• H1: Complicated grief symptom severity [tot_icg] will predict lower variability in spontaneous thought flow over time, as indexed by fewer functional connectivity state transitions over the resting state scan [statestats_numtransitions].

• H2: Complicated grief symptom severity [tot_icg] will predict greater automatic constraints on thought content, as indexed by more dwell time in states of default mode-salience network interconnectivity [statestats_dwelltime_s2].

• H3: Ruminative content in self-reported pre-resting state thought content [state_rum] will mediate the relationship between grief symptom severity [tot_icg] and automatic constraints as indexed by dwell time in states of default mode-salience network interconnectivity [statestats_dwelltime_s2].

Analysis
From the grant:

• H1 & H2: Two separate hierarchical regression models will be used to test (A) whether the continuous score of grief severity predicts dynamic functional connectivity variables (n transitions in H1; dwell time in H2), and (B) whether grief severity remains predictive when other plausible explanatory variables including but not limited to time since death, depressive symptom severity, and age are included in the model. Using R, I will manually build each linear regression model then perform model comparison to identify whether the model with the additional predictors performs significantly better than the simpler model, based on the F test.

For the power analysis, we had discussed just using the first model without the covariates, I think.

• H3: Observer-rated ruminative content in the post-scan reports…will be tested to identify whether the effect of grief severity on DFNC is explained by how much people were ruminating immediately before the resting state scan. A mediation model will be used to test the indirect and total effects of ruminative content on dwell time, using the R package “mediation”79.

# install.packages("mediation")
# install.packages("dplyr")
library(mediation)
library(dplyr)
filter <- dplyr::filter
dataforjf <- read_csv("/Users/sarenseeley/Dropbox/Grants\ and\ Applications/NRSA\ F31\ 2018-2019/Resubmission\ Reviews/power-analysis/dataforjf.csv")

Parsed with column specification:
cols(
  .default = col_double(),
  ID = col_character(),
  group = col_character(),
  tot_icg = col_integer(),
  tot_bdi = col_integer(),
  sex_m = col_character(),
  timesincedeath = col_integer(),
  tx_v1 = col_integer(),
  tx = col_integer(),
  meta_states_numstates = col_integer(),
  meta_states_totaldist = col_integer(),
  meta_states_changestates = col_integer(),
  meta_states_statespan = col_integer(),
  statestats_numtransitions = col_integer()
)
See spec(...) for full column specifications.

# Use PLACEBO SESSIONS ONLY for the following analyses
data_pl <- dataforjf %>% filter(tx == 0)
### Aim 1 Hypothesis 1
# First model with grief severity only
m1 <- lm(statestats_numtransitions ~ tot_icg, data = data_pl)
summary(m1)

Call:
lm(formula = statestats_numtransitions ~ tot_icg, data = data_pl)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3907 -2.2391  0.5366  1.8778  7.9650 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.71051    1.05108   6.384 2.14e-07 ***
tot_icg      0.03580    0.04053   0.883    0.383    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.091 on 36 degrees of freedom
Multiple R-squared:  0.02121,   Adjusted R-squared:  -0.005982 
F-statistic:  0.78 on 1 and 36 DF,  p-value: 0.383

# Second model including other hypothesized predictors 
# (here, depression symptoms + age + time since spouse's death at time of study):
m2 <- lm(statestats_numtransitions ~ tot_icg + tot_bdi + age_yrs + timesincedeath, data = data_pl) 
summary(m2)

Call:
lm(formula = statestats_numtransitions ~ tot_icg + tot_bdi + 
    age_yrs + timesincedeath, data = data_pl)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1019 -2.0838  0.0861  1.3641  6.9568 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)  
(Intercept)    11.2152610  5.2925972   2.119   0.0417 *
tot_icg        -0.0277563  0.0575273  -0.482   0.6326  
tot_bdi         0.1637882  0.0987070   1.659   0.1065  
age_yrs        -0.0698957  0.0757697  -0.922   0.3630  
timesincedeath  0.0003477  0.0020626   0.169   0.8672  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.047 on 33 degrees of freedom
Multiple R-squared:  0.1283,    Adjusted R-squared:  0.02265 
F-statistic: 1.214 on 4 and 33 DF,  p-value: 0.3234

# Comparing sum of squares & corresponding F statistics and p-values for m1 and m2: 
anova(m1,m2)

Analysis of Variance Table

Model 1: statestats_numtransitions ~ tot_icg
Model 2: statestats_numtransitions ~ tot_icg + tot_bdi + age_yrs + timesincedeath
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     36 344.02                           
2     33 306.38  3    37.644 1.3515 0.2746

### Aim 1 Hypothesis 2
# First model with grief severity only
m3 <- lm(statestats_dwelltime_s2 ~ tot_icg, data = data_pl)
summary(m3)

Call:
lm(formula = statestats_dwelltime_s2 ~ tot_icg, data = data_pl)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.095  -7.043  -2.783   3.531  30.975 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   0.6992     3.4777   0.201  0.84180   
tot_icg       0.3930     0.1341   2.930  0.00585 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.23 on 36 degrees of freedom
Multiple R-squared:  0.1926,    Adjusted R-squared:  0.1701 
F-statistic: 8.586 on 1 and 36 DF,  p-value: 0.00585

# Second model including other hypothesized predictors 
# (here, depression symptoms + age + time since spouse's death at time of study):
m4 <- lm(statestats_dwelltime_s2 ~ tot_icg + tot_bdi + age_yrs + timesincedeath, data = data_pl) 
summary(m4)

Call:
lm(formula = statestats_dwelltime_s2 ~ tot_icg + tot_bdi + age_yrs + 
    timesincedeath, data = data_pl)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.460  -6.495  -1.384   3.852  28.209 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)   
(Intercept)    -16.360323  17.465160  -0.937  0.35570   
tot_icg          0.615572   0.189835   3.243  0.00271 **
tot_bdi         -0.558718   0.325725  -1.715  0.09567 . 
age_yrs          0.205128   0.250034   0.820  0.41788   
timesincedeath   0.007230   0.006806   1.062  0.29584   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.05 on 33 degrees of freedom
Multiple R-squared:  0.2847,    Adjusted R-squared:  0.198 
F-statistic: 3.284 on 4 and 33 DF,  p-value: 0.02262

# Comparing sum of squares & corresponding F statistics and p-values for m3 and m4: 
anova(m3,m4)

Analysis of Variance Table

Model 1: statestats_dwelltime_s2 ~ tot_icg
Model 2: statestats_dwelltime_s2 ~ tot_icg + tot_bdi + age_yrs + timesincedeath
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     36 3766.2                           
2     33 3336.3  3    429.88 1.4174 0.2552

### Aim 1 Hypothesis 3
# model.M <- lm(M ~ X, myData)
# model.Y <- lm(Y ~ X + M, myData)
# results <- mediate(model.M, model.Y, treat='X', mediator='M',
#                   boot=TRUE, sims=100)
# X = tot_icg
# M = focus_spouse_perc
# Y = statestats_dwelltime_s2
mM <- lm(focus_spouse_perc ~ tot_icg, data_pl)
summary(mM)

Call:
lm(formula = focus_spouse_perc ~ tot_icg, data = data_pl)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46451 -0.15519 -0.00541  0.13085  0.64133 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 0.137369   0.088959   1.544  0.13129   
tot_icg     0.009622   0.003431   2.805  0.00807 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2616 on 36 degrees of freedom
Multiple R-squared:  0.1793,    Adjusted R-squared:  0.1565 
F-statistic: 7.866 on 1 and 36 DF,  p-value: 0.00807

mY <- lm(statestats_dwelltime_s2 ~ tot_icg + focus_spouse_perc, data_pl)
summary(mY)

Call:
lm(formula = statestats_dwelltime_s2 ~ tot_icg + focus_spouse_perc, 
    data = data_pl)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.932  -6.955  -2.662   4.076  30.950 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)   
(Intercept)         0.9679     3.6374   0.266  0.79172   
tot_icg             0.4118     0.1500   2.746  0.00946 **
focus_spouse_perc  -1.9564     6.5997  -0.296  0.76865   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.36 on 35 degrees of freedom
Multiple R-squared:  0.1946,    Adjusted R-squared:  0.1486 
F-statistic: 4.228 on 2 and 35 DF,  p-value: 0.02266

results <- mediation::mediate(mM, mY, treat='tot_icg', mediator='focus_spouse_perc',
                   boot=TRUE, sims=1000)
summary(results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value
ACME            -0.0188      -0.1588       0.0559    0.76
ADE              0.4118       0.0765       0.7380    0.00
Total Effect     0.3930       0.0819       0.7334    0.01
Prop. Mediated  -0.0479      -0.6164       0.2482    0.77

Sample Size Used: 38 


Simulations: 1000 

# ACME stands for Average Causal Mediation Effects
# ADE stands for Average Direct Effects
# Total Effect is a sum of a mediation (indirect) effect and a direct effect

Aim 2

Aim 2. To identify how intranasal oxytocin alters dynamic functional connectivity in older adults, and whether effects of oxytocin are moderated by grief severity.

• H1: Under oxytocin [tx == 1], the sample as a whole will show increased time in salience network-dominant states [statestats_dwelltime_s2], relative to placebo [tx == 0].

• H2: Complicated grief symptom severity will moderate oxytocin effects on DFNC [statestats_dwelltime_s2].

Analyses
From the grant:

• H1: Dependent t-tests will be used to compare oxytocin and placebo conditions to test whether oxytocin (A) increases dwell time in salience-network dominant states and/or (B) reduces n state transitions in the sample.

Hmm, not sure why I proposed doing t tests here - I think I would get the same info from the main effect of tx (oxytocin/placebo) in the model below for H2 because the observations are nested within person. Am I missing something here or was I just not thinking when I wrote this?

• H2: A multivariate linear mixed effects model will be used to test grief severity as a potential moderator of the relationship between condition (oxytocin vs. placebo) and DNFC dependent variables. Broadly, the two- intercept model80 will predict dwell time and n transitions as a function of “indicator variable” (i.e., which DV), condition, and ICG score (nested within person to account for interdependence), using maximum likelihood as implemented in the R package nlme81. This type of model allows for testing differential effects on the two DVs (dwell time; n transitions) and reduces chance of experiment-wise Type-I error. A multivariate linear mixed effects model also estimates the correlation between the DVs, and allows them to have different residual variances. It is most powerful when the dependent variables are correlated (as we might expect them to be here since dwell time and number of transitions are likely inversely related to each other).

I would revise this analysis to just do a regular linear mixed effects model with one outcome variable, statestats_dwelltime_s2. Dwell time in a particular state is more specific in terms of what it is measuring than number of transitions overall. I found a tutorial on multivariate random coefficient models using nlme but it’s a bit over my head at this time.

I also noted that dwell time in state 2 and n transitions are actually not correlated, as I expected them to be (corr.test(data$statestats_dwelltime_s2, data$statestats_numtransitions) outputs r = .05, p = .70).

library(nlme)

Attaching package: ‘nlme’

The following object is masked from ‘package:dplyr’:

    collapse

### Aim 2 Hypothesis 1
# This hypothesis tests the main effect of `tx` in the model below. I think it makes sense to have just random slopes and not random intercepts here but correct me if I'm wrong...
### Aim 2 Hypothesis 2
# This hypothesis tests the treatment x grief severity interaction in the model below.
m5 <- lme(statestats_dwelltime_s2 ~ tx*tot_icg, data = data, random = ~ 1|ID, method="ML")
summary(m5)

Linear mixed-effects model fit by maximum likelihood
 Data: data 
       AIC      BIC    logLik
  598.7121 612.6965 -293.3561

Random effects:
 Formula: ~1 | ID
        (Intercept) Residual
StdDev:      7.8943 9.141035

Fixed effects: statestats_dwelltime_s2 ~ tx * tot_icg 
                Value Std.Error DF    t-value p-value
(Intercept)  0.699176  4.219227 36  0.1657119  0.8693
tx          11.363483  4.515934 36  2.5163087  0.0165
tot_icg      0.392973  0.162709 36  2.4151859  0.0209
tx:tot_icg  -0.232347  0.174151 36 -1.3341660  0.1905
 Correlation: 
           (Intr) tx     tot_cg
tx         -0.535              
tot_icg    -0.879  0.470       
tx:tot_icg  0.470 -0.879 -0.535

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.5035893 -0.5997914 -0.1110756  0.2787453  2.6633116 

Number of Observations: 76
Number of Groups: 38