<style type="text/css">
div.main-container {
  max-width: 100%;
  margin-left: auto;
  margin-right: auto;
}
</style>

1. JP_Data

-Note: Large chunks of code have been hidden to improve readibility. -NOTE: Did some reading, this data is a particular type of MLM: “Cross-classified”. This means our level 1 units (i.e., RT observations) are clustered within two grouping variables (1 = Subject, 2 = Condition). Need to do a little more research on this. ## 1. Load Packages Used in Data Analysis -Code hidden

2. Data Read-In

-Code hidden

3. Data Pre-processing

Creating Trial Number Column

length(df$Session[df$Subject == 56001 & df$Session == 1]) #160 rows ("trials") per part.

## [1] 160

dim(df) #21,440 total rows ("trials")

## [1] 21440     8

21440/160 #160 rows repeated 134 times

## [1] 134

df$Trial <- rep(seq(1:160), 134) #Code creates a sequence of 1-160, per participant per session, and stores this in the "Trial" variable

Probe accuracy, incorrect trials, improbable reaction times

-The code below when run illustrates that there are no administrations #With <50% accuracy. Nothing excluded here.

df %>% 
  group_by(Subject, Session, Cue) %>% 
  summarise(Percent_Accurate = mean(Accuracy)) %>% 
  filter(Percent_Accurate < 0.50)

## `summarise()` has grouped output by 'Subject', 'Session'. You can override
## using the `.groups` argument.

## # A tibble: 0 × 4
## # Groups:   Subject, Session [0]
## # … with 4 variables: Subject <dbl>, Session <fct>, Cue <fct>,
## #   Percent_Accurate <dbl>

-The code below indicates the overall # of trials and the counts/percents for 1) incorrect trials to be excluded, and 2) improbable reaction time trials (i.e., <200ms OR >1500ms) to be excluded. The code immediately following then excludes this data.

df %>% 
  group_by(Session) %>% 
  summarise(Overall_Trials = n(), #Overall trial counts
            Count_Inc_Trial = sum(ifelse(Accuracy == 0,1,0)), #Incorrect trial count
            Per_Inc_Trial = 1 - mean(Accuracy), #Percent of trials that are incorrect
            Count_Inc_RT = sum(ifelse(RT < 200 | RT > 1500, 1, 0)), 
            #Count of trials whose reaction times are < 200 or > 1500
            Per_Inc_RT = 1 - mean(ifelse(RT < 200 | RT > 1500, 0, 1)))

## # A tibble: 6 × 6
##   Session Overall_Trials Count_Inc_Trial Per_Inc_Trial Count_Inc_RT Per_Inc_RT
##   <fct>            <int>           <dbl>         <dbl>        <dbl>      <dbl>
## 1 1                 4160              84        0.0202           79     0.0190
## 2 2                 4160             161        0.0387           80     0.0192
## 3 3                 3680              60        0.0163           70     0.0190
## 4 4                 3680              82        0.0223           62     0.0168
## 5 5                 2880              71        0.0247           69     0.0240
## 6 6                 2880              74        0.0257           90     0.0312

-Filtering out improbable reaction times

df <- df %>% 
  filter(Accuracy == 1) %>%  #Exclude incorrect trials
  filter(RT > 200 & RT < 1500) #Exclude RT's <200ms or >1500ms

-Maximum & Minimum RT after exclusions

max(df$RT)

## [1] 1499

min(df$RT)

## [1] 207

-Code for Median Absolute Deviation and SD calculations for Reaction Times hidden

-MAD & SD Calculation Results

df %>% 
  group_by(Session) %>% 
  summarise(Ovr_Trl_Aftr_Exc = n(), #Overall trial counts after excluding for above-noted criteria
            Count_MAD3_Over = sum(ifelse(RT > MAD_Exclude_Score, 1,0)), #Count of trials that had RT > 3MAD above individual median
            Per_MAD3_Over = mean(ifelse(RT > MAD_Exclude_Score, 1, 0)),
            Count_SD3_Over = sum(ifelse(RT > SD_Exclude_Score, 1,0)), #Count of trials that had RT > 3MAD above individual median
            Per_SD3_Over = mean(ifelse(RT > SD_Exclude_Score, 1, 0)))

## # A tibble: 6 × 6
##   Session Ovr_Trl_Aftr_Exc Count_MAD3_Over Per_MAD3_Over Count_SD3_Over Per_SD…¹
##   <fct>              <int>           <dbl>         <dbl>          <dbl>    <dbl>
## 1 1                   4005             205        0.0512             79   0.0197
## 2 2                   3931             138        0.0351             54   0.0137
## 3 3                   3551             115        0.0324             47   0.0132
## 4 4                   3538             119        0.0336             50   0.0141
## 5 5                   2746             114        0.0415             54   0.0197
## 6 6                   2725             128        0.0470             52   0.0191
## # … with abbreviated variable name ¹Per_SD3_Over

-Excluding Folks based on SD for Reaction Times

df <- df %>% 
  filter(RT < SD_Exclude_Score) #Exclude folks whose RT are >3SD above their individual means

3. Descriptive Statistics

df %>% select(RT, Condition, Session, Cue) %>%  group_by(Condition, Session, Cue) %>% 
  summarise(across(everything(), list(mean = mean, sd = sd, min = min, max = max)))

## `summarise()` has grouped output by 'Condition', 'Session'. You can override
## using the `.groups` argument.

## # A tibble: 36 × 7
## # Groups:   Condition, Session [18]
##    Condition Session Cue    RT_mean RT_sd RT_min RT_max
##    <chr>     <fct>   <fct>    <dbl> <dbl>  <dbl>  <dbl>
##  1 A         1       opioid    623.  210.    339   1447
##  2 A         1       pain      635.  213.    317   1440
##  3 A         2       opioid    580.  144.    207   1295
##  4 A         2       pain      576.  151.    367   1356
##  5 A         3       opioid    557.  125.    207   1167
##  6 A         3       pain      562.  133.    312   1316
##  7 A         4       opioid    551.  140.    342   1348
##  8 A         4       pain      552.  144.    328   1241
##  9 A         5       opioid    536.  141.    301   1217
## 10 A         5       pain      534.  137.    322   1219
## # … with 26 more rows

4. Attentional Bias (AB) Calculations

-Code for AB calculations hidden

-Attentional Bias Metrics Per Subject, Condition, Session, & Cue

By_Subject_Output

## # A tibble: 268 × 10
## # Groups:   Subject, Condition, Session [134]
##    Subject Condi…¹ Session Cue   Mean_…² Mean_…³ Mean_…⁴ Peak_…⁵ Peak_…⁶ Varia…⁷
##      <dbl> <chr>   <fct>   <fct>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1   56001 C       1       opio…   6.24     97.9   103.      442     203    82.1
##  2   56001 C       1       pain   33.2     139.     61.2     429     226    76.5
##  3   56001 C       2       opio…  12.3     114.    110.      402     296    87.9
##  4   56001 C       2       pain    0.776   100.    107.      380     287    76.7
##  5   56001 A       3       opio…  22.9     115.    113.      350     528   101. 
##  6   56001 A       3       pain    4.97     65.3    89.6     228     228    74.2
##  7   56001 A       4       opio… -39.8      82.1   106.      240     347    73.7
##  8   56001 A       4       pain    1.86     75.7   101.      326     486    80.2
##  9   56001 B       5       opio…  -0.504    88.8    65.7     251     178    64.5
## 10   56001 B       5       pain    9.94    110.     52.1     305     289    79.9
## # … with 258 more rows, and abbreviated variable names ¹Condition, ²Mean_Bias,
## #   ³Mean_Toward, ⁴Mean_Away, ⁵Peak_Toward, ⁶Peak_Away, ⁷Variability

By_Condition_Output

## # A tibble: 36 × 9
## # Groups:   Condition, Session [18]
##    Condition Session Cue    Mean_Bias Mean_Tow…¹ Mean_…² Peak_…³ Peak_…⁴ Varia…⁵
##    <chr>     <fct>   <fct>      <dbl>      <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 A         1       opioid     -4.88      112.    117.     323.    317.    92.8
##  2 A         1       pain        1.90      117.    125.     400.    353.    98.4
##  3 B         1       opioid     -3.52       95.1    90.3    318.    313.    89.7
##  4 B         1       pain        6.39       84.6    90.5    292.    290.    81.8
##  5 C         1       opioid     -9.58      119.    128.     364.    367.   104. 
##  6 C         1       pain        5.10      130.    118.     396.    386.   107. 
##  7 A         2       opioid     20.0       102.     96.6    317.    329.    89.9
##  8 A         2       pain      -10.2        93.9   105.     354.    341.    89.9
##  9 B         2       opioid      3.00      106.    109.     406.    331.    86.2
## 10 B         2       pain       28.8       117.     96.1    364.    299.    94.3
## # … with 26 more rows, and abbreviated variable names ¹Mean_Toward, ²Mean_Away,
## #   ³Peak_Toward, ⁴Peak_Away, ⁵Variability

By_Session_Output

## # A tibble: 12 × 8
## # Groups:   Session [6]
##    Session Cue    Mean_Bias Mean_Toward Mean_Away Peak_Toward Peak_Away Variab…¹
##    <fct>   <fct>      <dbl>       <dbl>     <dbl>       <dbl>     <dbl>    <dbl>
##  1 1       opioid    -6.27         109.      113.        337.      335.     96.2
##  2 1       pain       4.50         113.      112.        366.      347.     96.9
##  3 2       opioid     1.82         108.      110.        357.      341.     93.2
##  4 2       pain       6.25         110.      112.        359.      342.     96.2
##  5 3       opioid    -0.233        111.      108.        373.      364.     94.1
##  6 3       pain      -4.52         105.      111.        338.      353.     94.5
##  7 4       opioid     4.76         111.      110.        377.      352.     93.3
##  8 4       pain       2.01         110.      106.        343.      362.     96.2
##  9 5       opioid    -3.19         106.      105.        324.      348.     92.1
## 10 5       pain     -13.3          103.      107.        326.      342.     90.4
## 11 6       opioid    -3.74         105.      110.        335.      329.     96.1
## 12 6       pain     -12.3          106.      116.        313.      348.     91.9
## # … with abbreviated variable name ¹Variability

By_Cue_Output

## # A tibble: 2 × 7
##   Cue    Mean_Bias Mean_Toward Mean_Away Peak_Toward Peak_Away Variability
##   <fct>      <dbl>       <dbl>     <dbl>       <dbl>     <dbl>       <dbl>
## 1 opioid     -1.01        108.      110.        352.      345.        94.2
## 2 pain       -1.76        108.      111.        343.      349.        94.7

5. ANOVA’s

-Code Hidden -Changing session variable to represent just pre-post

-Something to note: There are several subjects who appear to have no data for at least Condition/1 time-point:

Folks Missing 1 pre/post data point: #Condition A: 15 #Condition B: 02 #Condition C: 14

Folks Missing both pre/post data points #Condition A: 03, 11, 13, 14, 18, 22, 28 #Condition B: 03, 13, 17 #Condition C: 17,18,23

R seems to be dropping these folks from ANOVA analysis below. I’m trying to specify a way for them to be kept and not succeeding at the moment, but even if I can, that still is going to result with at least 13/28 getting dropped. A REALLY small sample of N=15 at best or N=12 at worst.

6a. Adding grouping by condition variable

df$Dose <- ifelse(df$Subject == "56001"  | df$Subject == "56002"| df$Subject == "56005"| df$Subject == "56006" | df$Subject == "56009"| df$Subject == "56010"| df$Subject == "56014"| df$Subject == "56016"| df$Subject == "56018"| df$Subject == "56019"| df$Subject == "56020"| df$Subject == "56021"| df$Subject == "56022"| df$Subject == "56023"| df$Subject == "56026"| df$Subject == "56027", "Low", "High")

df$Dose[df$Subject == "56024"] <- "Not randomized"
table(df$Dose)

## 
##  High   Low 
##  7607 12553

-Changing Data into Long Format

df_long <- df %>% select(Subject, Dose, Condition, Session, Cue, mean_bias:variability) %>% 
  pivot_longer(cols = c(mean_bias:variability), names_to = "Bias_Type", values_to = "RT") %>% 
  mutate(Subject = factor(Subject)) %>% 
  distinct()

dim(df_long) #Dataframe dimensions

## [1] 1608    7

# Number of RT's each participant contributes to the dataframe
table(df_long$Subject)

## 
## 56001 56002 56003 56004 56005 56006 56007 56008 56009 56010 56011 56012 56013 
##    72    60    36    72    72    72    72    72    72    72    48    72    24 
## 56014 56015 56016 56017 56018 56019 56020 56021 56022 56023 56025 56026 56027 
##    36    60    72    24    24    72    72    72    48    48    72    72    72 
## 56028 
##    48

data.frame(table(df_long$Subject)) %>% summarise(Min = min(Freq), Max = max(Freq))

##   Min Max
## 1  24  72

6. ANOVAs

6a. Mean Bias

-DV: Opioid Mean Bias

-DV: Pain Mean Bias

6b. Mean Away Bias

-DV: Opioid Mean Away Bias

-DV: Pain Mean Away Bias

6c. Mean Toward Bias

-DV: Opioid Mean Toward Bias

-DV: Pain Mean Toward Bias

6d. Peak Away Bias

-DV: Opioid Peak Away Bias

-DV: Pain Peak Away Bias

6e. Peak Toward Bias

-DV: Opioid Peak Toward Bias

-DV: Pain Peak Toward Bias

6f. Variability Bias

-DV: Opioid Variability Bias

-DV: Pain Variability Bias

7. Muiltilevel Models

7a. Mean Bias

-DV: Opioid Mean Bias

re_O_MB_LD <- lmer(RT ~ (1|Subject), data = df_long_opioid_mean_bias %>% filter(Dose == "Low"))
re_O_MB_HD <- lmer(RT ~ (1|Subject), data = df_long_opioid_mean_bias %>% filter(Dose == "High"))

performance::icc(re_O_MB_LD)

## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.043
##   Conditional ICC: 0.043

performance::icc(re_O_MB_HD)

## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.231
##   Conditional ICC: 0.231

MLM_O_MB_LD <- lmer(RT ~ Session*Condition + (1 + Session|Subject), data = df_long_opioid_mean_bias %>% filter(Dose == "Low"))
MLM_O_MB_HD <- lmer(RT ~ Session*Condition + (1 + Session|Subject), data = df_long_opioid_mean_bias %>% filter(Dose == "High"))

## boundary (singular) fit: see help('isSingular')

#The code below then obtains mltiple Degree-of-Freedom F-Tests for the above-specified model
anova(MLM_O_MB_LD, type = "III")

## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## Session            187.78  187.78     1 15.885  0.2607 0.6166
## Condition          199.17   99.58     2 51.463  0.1383 0.8712
## Session:Condition 1545.79  772.89     2 50.360  1.0731 0.3496

anova(MLM_O_MB_HD, type = "III")

## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## Session            193.58  193.58     1 11.849  0.3312 0.5757
## Condition         2612.92 1306.46     2 38.010  2.2353 0.1208
## Session:Condition  156.79   78.39     2 35.556  0.1341 0.8749

#Model Parameters 
summary(MLM_O_MB_LD)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ Session * Condition + (1 + Session | Subject)
##    Data: df_long_opioid_mean_bias %>% filter(Dose == "Low")
## 
## REML criterion at convergence: 769.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8170 -0.6311 -0.0437  0.4117  3.4516 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  Subject  (Intercept) 507.9    22.54         
##           SessionPost 840.4    28.99    -0.96
##  Residual             720.2    26.84         
## Number of obs: 84, groups:  Subject, 16
## 
## Fixed effects:
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)              -5.894      9.514  39.086  -0.620    0.539
## SessionPost               8.037     13.003  41.673   0.618    0.540
## ConditionB               10.530     10.402  51.274   1.012    0.316
## ConditionC               -2.869     10.504  50.044  -0.273    0.786
## SessionPost:ConditionB  -14.898     14.530  50.744  -1.025    0.310
## SessionPost:ConditionC    5.206     14.929  49.672   0.349    0.729
## 
## Correlation of Fixed Effects:
##             (Intr) SssnPs CndtnB CndtnC SsP:CB
## SessionPost -0.798                            
## ConditionB  -0.597  0.442                     
## ConditionC  -0.578  0.426  0.529              
## SssnPst:CnB  0.432 -0.619 -0.721 -0.382       
## SssnPst:CnC  0.408 -0.589 -0.374 -0.707  0.527

summary(MLM_O_MB_HD)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ Session * Condition + (1 + Session | Subject)
##    Data: df_long_opioid_mean_bias %>% filter(Dose == "High")
## 
## REML criterion at convergence: 426
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.18606 -0.50718 -0.09435  0.55229  3.02205 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  Subject  (Intercept)  53.0     7.28        
##           SessionPost 103.1    10.15    1.00
##  Residual             584.5    24.18        
## Number of obs: 50, groups:  Subject, 11
## 
## Fixed effects:
##                        Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)               2.210     10.208  30.630   0.217    0.830
## SessionPost               9.670     13.911  30.831   0.695    0.492
## ConditionB                3.511     12.844  35.533   0.273    0.786
## ConditionC              -12.442     12.600  35.525  -0.987    0.330
## SessionPost:ConditionB   -9.120     17.963  35.504  -0.508    0.615
## SessionPost:ConditionC   -6.626     17.387  35.396  -0.381    0.705
## 
## Correlation of Fixed Effects:
##             (Intr) SssnPs CndtnB CndtnC SsP:CB
## SessionPost -0.628                            
## ConditionB  -0.758  0.538                     
## ConditionC  -0.775  0.545  0.618              
## SssnPst:CnB  0.527 -0.731 -0.698 -0.425       
## SssnPst:CnC  0.538 -0.763 -0.426 -0.699  0.587
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

#Means/SDs
df_long_opioid_mean_bias %>% filter(Dose == "Low") %>% select(-Cue, -Bias_Type) %>% group_by(Session, Condition) %>% summarise(M = mean(RT), SD = sd(RT))

## `summarise()` has grouped output by 'Session'. You can override using the
## `.groups` argument.

## # A tibble: 6 × 4
## # Groups:   Session [2]
##   Session Condition      M    SD
##   <fct>   <chr>      <dbl> <dbl>
## 1 Pre     A          -6.36  35.2
## 2 Pre     B           4.70  30.8
## 3 Pre     C         -10.6   32.4
## 4 Post    A           2.36  37.3
## 5 Post    B          -2.23  21.0
## 6 Post    C           4.81  35.5

df_long_opioid_mean_bias %>% filter(Dose == "High") %>% select(-Cue, -Bias_Type) %>% group_by(Session, Condition) %>% summarise(M = mean(RT), SD = sd(RT))

## `summarise()` has grouped output by 'Session'. You can override using the
## `.groups` argument.

## # A tibble: 6 × 4
## # Groups:   Session [2]
##   Session Condition       M    SD
##   <fct>   <chr>       <dbl> <dbl>
## 1 Pre     A           0.291  23.2
## 2 Pre     B           5.44   17.5
## 3 Pre     C         -10.7    35.9
## 4 Post    A          11.5    30.1
## 5 Post    B           5.75   29.1
## 6 Post    C          -8.41   25.2

emmeans::emmeans(MLM_O_MB_LD, specs = c("Session", "Condition"))

##  Session Condition emmean   SE   df lower.CL upper.CL
##  Pre     A          -5.89 9.56 38.7    -25.2     13.4
##  Post    A           2.14 7.96 52.4    -13.8     18.1
##  Pre     B           4.64 8.99 34.8    -13.6     22.9
##  Post    B          -2.23 7.13 50.2    -16.5     12.1
##  Pre     C          -8.76 9.27 36.5    -27.6     10.0
##  Post    C           4.48 7.96 52.4    -11.5     20.5
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

emmeans::emmeans(MLM_O_MB_HD, specs = c("Session", "Condition"))

##  Session Condition emmean    SE   df lower.CL upper.CL
##  Pre     A           2.21 10.72 31.2    -19.6    24.06
##  Post    A          11.88 11.19 27.1    -11.1    34.83
##  Pre     B           5.72  8.51 29.6    -11.7    23.11
##  Post    B           6.27 10.52 25.1    -15.4    27.94
##  Pre     C         -10.23  8.02 29.1    -26.6     6.18
##  Post    C          -7.19  9.42 22.1    -26.7    12.35
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

#Actual Means
AM_LD <- ggplot(df_long_opioid_mean_bias %>% filter(Dose == "Low"), aes(x=Session,y=RT, col = Condition)) +
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5)) + ggtitle("Group Means") +
  theme_light() +
  theme(plot.title = element_text(hjust = 0.5))
AM_HD <- ggplot(df_long_opioid_mean_bias %>% filter(Dose == "High"), aes(x=Session,y=RT, col = Condition)) +
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5)) + ggtitle("Group Means") +
  theme_light() +
  theme(plot.title = element_text(hjust = 0.5))

#Predicted Means
PM_LD <- ggplot(df_long_opioid_mean_bias %>% filter(Dose == "Low"), aes(x=Session,y=predict(MLM_O_MB_LD), col = Condition)) +
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5)) + ggtitle("Group Means") +
  theme_light() +
  theme(plot.title = element_text(hjust = 0.5))
PM_HD <- ggplot(df_long_opioid_mean_bias %>% filter(Dose == "High"), aes(x=Session,y=predict(MLM_O_MB_HD), col = Condition)) +
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5)) + ggtitle("Group Means") +
  theme_light() +
  theme(plot.title = element_text(hjust = 0.5))

cowplot::plot_grid(AM_LD,AM_HD)

cowplot::plot_grid(PM_LD,PM_HD)