SOCIAL CONNECTION SCALE ’

Items designed to assess social connection were collected for 5 days in a row via participants smartphones at 1-week post intervention and 14-weeks post intervention.

Two composite scores were created: 1) Positive, close, and connected social interactions (Total number of social interactions, Number different people interacted with, Social closeness/connection, Number of positive social interactions, and Degree of positive social interactions. 2) Negative Interactions (Number of negative social interactions and Degree of negative social interactions)

Exclusion Criteria (per scale):

Step 1: If Pp has complete or partial data for 15% of EMA timepoints > included

Step 2: If Pp has completed 15% of all item responses > included

Step 4: Series means will be imputed for any remaining missing data

Step 5: Create composite score with included, imputed items

Evaluating missingness for state social connection items

Step 1 - Percentage of participants who are excluded based on total timepoints completed

library(dplyr)

# Extract column for the daily social connection scale
socialconnect_scale <- eq_EMA_data %>%
  dplyr::select(contains("soc_"))

socialconnect_scale$PID <- eq_EMA_data$PID
socialconnect_scale$EMA_timepoint <- eq_EMA_data$EMA_timepoint
socialconnect_scale$Condition <- eq_EMA_data$Condition

#socialconnect_scale

#Count the total number of timepoints completed per participant
socialconnect_NumComplete <- count(socialconnect_scale, vars = PID)
socialconnect_NumComplete$PercentComplete <- socialconnect_NumComplete$n / 14
colnames(socialconnect_NumComplete)[1:2] <- c("PID", "surveyComplete")

# Extract the PID with less than 15% completion rate across all timepoints
socialconnect_exclude1 <- socialconnect_NumComplete$PID[socialconnect_NumComplete$PercentComplete < 0.15]

#socialconnect_exclude1

# Calculate the percentage of participants who are excluded based on total timepoints completed
socialconnect_percent_AllTime <- length(socialconnect_exclude1)/113
socialconnect_percent_AllTime
## [1] 0.01769912
#Remove participants from the scale dataset who have been excluded based on the first criteria
socialconnect_scale_include1 <- socialconnect_scale %>%
  filter(! PID %in% socialconnect_exclude1)

Step 2 - Percentage of remaining participants excluded due to less than 15% of all item responses completed

# Calculate the number of total completed data per row
socialconnect_scale_include1$itemComplete <- rowSums(!is.na(socialconnect_scale_include1[,1:7]))

#socialconnect_scale_include1

# Aggregate within subject
socialconnect_itemComplete <- aggregate(itemComplete ~ PID, data = socialconnect_scale_include1 [,8:11], sum)

#socialconnect_itemComplete

# Calculate the percentage of missing
socialconnect_totalItem <- 7 * 14 #4 items * 14 timepoints 
socialconnect_itemComplete$percentComplete <- round(socialconnect_itemComplete$itemComplete/socialconnect_totalItem,3) #round to 3 decimal points

#socialconnect_itemComplete #NEEDS UPDATED - MORE THAN 100% COMPLETE - WHY?

# Extract the PID with less than 15% completion rate across both row and column 
#Brendan suggests higher threshold here bc it is related to the reliability of the items. You can systematically look at the reliability and plot it as a functions of items. 
socialconnect_exclude2 <- socialconnect_itemComplete$PID[socialconnect_itemComplete$percentComplete < 0.15]

#socialconnect_exclude2 

# Calculate the percentage of participant who are excluded based on this scale
length(socialconnect_exclude2)/113
## [1] 0
#Remove participants from the scale dataset who have been excluded based on the second criteria
socialconnect_scale_include2 <- socialconnect_scale_include1 %>%
  filter(! PID %in% socialconnect_exclude2)

#socialconnect_scale_include2

Step 3: If a scale has 70% of participant responses after previous steps, scale included, Count the percentage of participants remaining after step 1 and 2 who are included in the scale.

# Count the percentage of participants remaining after step 1 and 2 who are included in the scale
totalIncluded <- length(unique(socialconnect_scale_include2$PID))
percentageIncluded <- round((totalIncluded/113),3)
percentageIncluded
## [1] 0.982

Step 4: Impute series means for any remaining missing data

socialconnect_scale_include2$soc_numdifpeople_imp <- ifelse(is.na(socialconnect_scale_include2$soc_numdifpeople), mean(socialconnect_scale_include2$soc_numdifpeople, na.rm=TRUE), socialconnect_scale_include2$soc_numdifpeople)

socialconnect_scale_include2$soc_numpeople_imp <- ifelse(is.na(socialconnect_scale_include2$soc_numpeople), mean(socialconnect_scale_include2$soc_numpeople, na.rm=TRUE), socialconnect_scale_include2$soc_numpeople)

socialconnect_scale_include2$soc_positive_imp <- ifelse(is.na(socialconnect_scale_include2$soc_positive), mean(socialconnect_scale_include2$soc_positive, na.rm=TRUE), socialconnect_scale_include2$soc_positive)

socialconnect_scale_include2$soc_neg_imp <- ifelse(is.na(socialconnect_scale_include2$soc_neg), mean(socialconnect_scale_include2$soc_neg, na.rm=TRUE), socialconnect_scale_include2$soc_neg)

socialconnect_scale_include2$soc_howpos_imp <- ifelse(is.na(socialconnect_scale_include2$soc_howpos), 
mean(socialconnect_scale_include2$soc_howpos, na.rm=TRUE), 
socialconnect_scale_include2$soc_howpos) #WHY ISN'T THIS ONE IMPUTING?

socialconnect_scale_include2$soc_howneg_imp <- ifelse(is.na(socialconnect_scale_include2$soc_howneg), mean(socialconnect_scale_include2$soc_howneg, na.rm=TRUE), socialconnect_scale_include2$soc_howneg)

socialconnect_scale_include2$soc_closeness_imp <- ifelse(is.na(socialconnect_scale_include2$soc_closeness), mean(socialconnect_scale_include2$soc_closeness, na.rm=TRUE), socialconnect_scale_include2$soc_closeness)

Step 5: Creating composite scores for Social Connection

Re-coding score “soc_closeness” - score will now reflect that higher scores on soc_closeness_imp will equate to greater social connectedness.

socialconnect_scale_include2$soc_closeness_imp <- 5-socialconnect_scale_include2$soc_closeness_imp

Correlation Matrix - All items

socialconnect_scale_include_cor <- dplyr::select(socialconnect_scale_include2, soc_numdifpeople_imp:soc_closeness_imp)

# All Items
sapply(socialconnect_scale_include_cor, is.factor)
## soc_numdifpeople_imp    soc_numpeople_imp     soc_positive_imp 
##                FALSE                FALSE                FALSE 
##          soc_neg_imp       soc_howpos_imp       soc_howneg_imp 
##                FALSE                FALSE                FALSE 
##    soc_closeness_imp 
##                FALSE
cor <- round(cor(socialconnect_scale_include_cor[sapply(socialconnect_scale_include_cor, function(x) !is.factor(x))], use = "complete.obs"),2) #Correlation matrix of indices 

#cor

library(reshape2)
melted_cormat <- melt(cor)
#head(melted_cormat)

# Get lower triangle of the correlation matrix
  get_lower_tri<-function(cormat){
    cormat[upper.tri(cormat)] <- NA
    return(cormat)
  }
  # Get upper triangle of the correlation matrix
  get_upper_tri <- function(cormat){
    cormat[lower.tri(cormat)]<- NA
    return(cormat)
  }
  
  upper_tri <- get_upper_tri(cor)
#upper_tri
  
  # Melt the correlation matrix
library(reshape2)
melted_cormat <- melt(upper_tri, na.rm = TRUE)


reorder_cormat <- function(cor){
# Use correlation between variables as distance
dd <- as.dist((1-cor)/2)
hc <- hclust(dd)
cor <-cor[hc$order, hc$order]
}

# Reorder the correlation matrix
cor <- reorder_cormat(cor)
upper_tri <- get_upper_tri(cor)
# Melt the correlation matrix
melted_cormat <- melt(upper_tri, na.rm = TRUE)
# Create a ggheatmap
library(ggplot2)
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
 geom_tile(color = "white")+
 scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
   midpoint = 0, limit = c(-1,1), space = "Lab", 
    name="Pearson\nCorrelation") +
  theme_minimal()+ # minimal theme
 theme(axis.text.x = element_text(angle = 45, vjust = 1, 
    size = 12, hjust = 1))+
 coord_fixed()


cor_heatmap <- ggheatmap + 
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme(
  axis.title.x = element_blank(),
  axis.title.y = element_blank(),
  panel.grid.major = element_blank(),
  panel.border = element_blank(),
  panel.background = element_blank(),
  axis.ticks = element_blank(),
  legend.justification = c(1, 0),
  legend.position = c(0.6, 0.7),
  legend.direction = "horizontal")+
  guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
                title.position = "top", title.hjust = 0.5))

cor_heatmap

Looking at the correlation matrix for the social connection items, we see that soc_neg_imp and soc_howneg_imp do not correlate well with the other scale items (but do correlate well with one another). We will create a composite of the remaining items and then a composite of the negative social connection items. (see below).

Items / Composites that will be analyzed

Positive, close, and connected social interactions

1a - Total number of social interactions (soc_numpeople_imp) (N/N/N)

1b - Number different people interacted with (soc_numdifpeople_imp) (N/Y/Y)

2 - Social closeness/connection (soc_closeness_imp) (N/Y/Y)(N/Y/Y)

3a - Number of positive social interactions (soc_positive_imp) (N/N/N)

3b - Degree of positive social interactions (soc_howpos_imp) (N/N/N)

Negative Interactions

3c - Number of negative social interactions (soc_neg_imp) (N/N/Y)

3d - Degree of negative social interactions (soc_howneg_imp) (N/N/N)

4 - Composite of positive, close, and connected social interactions items (EMA_social_comp) (N/Y/Y)

5 - Composite of negative interactions (EMA_social_comp_neg) (N/Y/Y)

Checking for internal validity of Social connection composite and creating composite score.

#Checking for internal validity of positive, close, and connected social connection composite 

df_social_connect_comp <- dplyr::select(socialconnect_scale_include2, soc_numdifpeople_imp, soc_numpeople_imp, soc_positive_imp, soc_howpos_imp, soc_closeness_imp) 
                            
psych::alpha(df_social_connect_comp) #calculating chronbach's alpha for pos social connect items
## 
## Reliability analysis   
## Call: psych::alpha(x = df_social_connect_comp)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N  ase mean   sd median_r
##       0.74      0.79     0.8      0.43 3.7 0.01  2.5 0.91     0.38
## 
##  lower alpha upper     95% confidence boundaries
## 0.72 0.74 0.76 
## 
##  Reliability if an item is dropped:
##                      raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## soc_numdifpeople_imp      0.67      0.73    0.74      0.40 2.7   0.0126 0.036
## soc_numpeople_imp         0.65      0.72    0.70      0.39 2.5   0.0128 0.023
## soc_positive_imp          0.62      0.69    0.70      0.36 2.2   0.0149 0.042
## soc_howpos_imp            0.75      0.82    0.81      0.52 4.4   0.0115 0.041
## soc_closeness_imp         0.80      0.78    0.78      0.46 3.5   0.0067 0.070
##                      med.r
## soc_numdifpeople_imp  0.38
## soc_numpeople_imp     0.38
## soc_positive_imp      0.30
## soc_howpos_imp        0.54
## soc_closeness_imp     0.47
## 
##  Item statistics 
##                         n raw.r std.r r.cor r.drop mean   sd
## soc_numdifpeople_imp 1345  0.75  0.78  0.73   0.62  3.0 1.03
## soc_numpeople_imp    1345  0.78  0.80  0.79   0.64  3.2 1.19
## soc_positive_imp     1345  0.84  0.85  0.83   0.73  2.6 1.24
## soc_howpos_imp       1345  0.49  0.58  0.41   0.39  2.2 0.53
## soc_closeness_imp    1345  0.76  0.68  0.54   0.44  1.7 2.01
alpha_connect <- psych::alpha(df_social_connect_comp)$total$std.alpha #extracting just the alpha value for items  

alpha_connect <- round(alpha_connect, 2)  

print(alpha_connect) #Good internal consistency
## [1] 0.79
#Creating column with social connection composite score ("EMA_social_comp")

socialconnect_scale_include2$EMA_social_comp <-socialconnect_scale_include2$soc_numdifpeople_imp + socialconnect_scale_include2$soc_numpeople_imp + socialconnect_scale_include2$soc_positive_imp + socialconnect_scale_include2$soc_howpos_imp + socialconnect_scale_include2$soc_closeness_imp

Write up- Internal validitiy of social connection

Checking for internal validity of positive, close, and connected social connection composite reveals a chronbach’s alpha (a= .79).

Checking for internal validity of positive social interaction composite.

#Checking for internal validity of positive, close, and connected social connection composite without soc_numdifpeople_imp, soc_numpeople_imp 

df_social_positive <- dplyr::select(socialconnect_scale_include2, soc_howpos_imp, soc_closeness_imp) 
                            
psych::alpha(df_social_positive) #calculating chronbach's alpha for pos social connect items
## 
## Reliability analysis   
## Call: psych::alpha(x = df_social_positive)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean  sd median_r
##       0.36      0.61    0.44      0.44 1.6 0.016    2 1.1     0.44
## 
##  lower alpha upper     95% confidence boundaries
## 0.32 0.36 0.39 
## 
##  Reliability if an item is dropped:
##                   raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## soc_howpos_imp         0.44      0.44     0.2      0.44  NA       NA  0.44
## soc_closeness_imp      0.20      0.44      NA        NA  NA       NA  0.20
##                   med.r
## soc_howpos_imp     0.44
## soc_closeness_imp  0.44
## 
##  Item statistics 
##                      n raw.r std.r r.cor r.drop mean   sd
## soc_howpos_imp    1345  0.62  0.85  0.56   0.44  2.2 0.53
## soc_closeness_imp 1345  0.98  0.85  0.56   0.44  1.7 2.01
## 
## Non missing response frequency for each item
##                     -3   -2   -1    1 1.73033707865169    2 2.21062992125984
## soc_howpos_imp    0.00 0.00 0.00 0.08             0.00 0.45             0.24
## soc_closeness_imp 0.07 0.02 0.09 0.28             0.01 0.00             0.00
##                      3    4 miss
## soc_howpos_imp    0.23 0.00    0
## soc_closeness_imp 0.37 0.17    0
alpha_positive <- psych::alpha(df_social_positive)$total$std.alpha #extracting just the alpha value for items  

alpha_positive <- round(alpha_positive, 2)  

print(alpha_positive) #Good internal consistency
## [1] 0.61
#Creating column with positive social connection composite score ("EMA_social_comp")

socialconnect_scale_include2$EMA_social_pos <- socialconnect_scale_include2$soc_positive_imp + socialconnect_scale_include2$soc_howpos_imp + socialconnect_scale_include2$soc_closeness_imp

Write up- Internal validity of positive social connection

Checking for internal validity of positive social connection composite reveals a chronbach’s alpha (a= .61).

Checking for internal validity of negative social interaction composite.

#Checking for internal validity of Negative social connection composite 

df_neg_social_connect_comp <- dplyr::select(socialconnect_scale_include2, soc_neg_imp, soc_howneg_imp) 
                            
psych::alpha(df_neg_social_connect_comp) #calculating chronbach's alpha for neg social connect items
## 
## Reliability analysis   
## Call: psych::alpha(x = df_neg_social_connect_comp)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
##       0.55       0.6    0.43      0.43 1.5 0.021 0.93 0.58     0.43
## 
##  lower alpha upper     95% confidence boundaries
## 0.51 0.55 0.59 
## 
##  Reliability if an item is dropped:
##                raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## soc_neg_imp         0.43      0.43    0.19      0.43  NA       NA  0.43  0.43
## soc_howneg_imp      0.19      0.43      NA        NA  NA       NA  0.19  0.43
## 
##  Item statistics 
##                   n raw.r std.r r.cor r.drop mean   sd
## soc_neg_imp    1345  0.92  0.85  0.56   0.43  0.6 0.84
## soc_howneg_imp 1345  0.75  0.85  0.56   0.43  1.3 0.51
## 
## Non missing response frequency for each item
##                   0 0.602400600150037    1 1.25987841945289    2    3    4 5
## soc_neg_imp    0.57              0.01 0.30             0.00 0.09 0.03 0.01 0
## soc_howneg_imp 0.00              0.00 0.76             0.02 0.18 0.04 0.00 0
##                miss
## soc_neg_imp       0
## soc_howneg_imp    0
alpha_neg <- psych::alpha(df_neg_social_connect_comp)$total$std.alpha #extracting just the alpha value for items  

alpha_neg <- round(alpha_neg, 2)  

print(alpha_neg) #Acceptable internal consistency
## [1] 0.6
#Creating column with social connection composite score ("EMA_social_comp")

socialconnect_scale_include2$EMA_neg_social_comp <-socialconnect_scale_include2$soc_neg_imp + socialconnect_scale_include2$soc_howneg_imp

Write up- Internal validity of negative social connection

Checking for internal validity of negative social connection composite reveals a chronbach’s alpha (a= .60)

Re-labeling condition labels. The conditions were originally labeled as 1 and 2 in the dataset, this function changes this such that condition 1 is now condition 0 and condition 2 is now condition 1.

library(tidyr)
socialconnect_scale_include2 <- socialconnect_scale_include2 %>% 
  separate (EMA_timepoint, into = c("onewkpost_14wkpost", "timepoint")) %>% 
  mutate(Condition = Condition -1) #the conditions were originally labeled as 1 and 2 in the dataset, this function changes this such that condition 1 is now condition 0 and condition 2 is now condition 1. 
#Note: 
  # Original condition 1 in the dataset is now "Condition 0"
  # Original condition 2 in the dataset is now "Condition 1" #this can be written more clearly with an ifelse statement

#socialconnect_scale_include2

Creating separate dataframes for 1-week and 14-week follow-up

socialconnect_scale_1wkpost <- socialconnect_scale_include2 %>% 
  filter(onewkpost_14wkpost == '1week')

socialconnect_scale_14wkpost <- socialconnect_scale_include2 %>% 
  filter(onewkpost_14wkpost == '14week')

Re-shaping variables to factors

library(lme4)
library(lmerTest)

socialconnect_scale_include2 <- socialconnect_scale_include2 %>% 
  mutate(onewkpost_14wkpost = as.factor(onewkpost_14wkpost), 
         timepoint = as.factor(timepoint))

Mixed models for EMA Social Connection Scales

POSITIVE, CLOSE, and CONNECTED SOCIAL CONNECTION COMPOSITE

Describing the data

# Descriptives - EMA Social Connection Composite (EMA_social_comp)

#BY TIME
# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$EMA_social_comp)
##    vars   n  mean   sd median trimmed  mad min max range  skew kurtosis   se
## X1    1 670 13.23 4.16  13.61   13.49 3.86   0  22    22 -0.54    -0.01 0.16
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$EMA_social_comp)
##    vars   n  mean   sd median trimmed  mad min max range  skew kurtosis   se
## X1    1 675 12.13 4.83  12.21    12.3 4.76  -2  22    24 -0.37     0.01 0.19
#BY CONDITION / TIME
# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$EMA_social_comp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 323 13.57223 3.910311     14 13.80629 4.4478   0  22    22
## X12    2      1    1 347 12.90287 4.359070     13 13.17271 4.4478   1  22    21
##           skew    kurtosis        se
## X11 -0.5183477  0.07009815 0.2175755
## X12 -0.5203587 -0.16305883 0.2340071
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$EMA_social_comp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 344 12.65017 4.463040     13 12.74781 4.4478  -2  22    24
## X12    2      1    1 331 11.58084 5.138869     12 11.79064 5.9304  -2  22    24
##           skew   kurtosis        se
## X11 -0.2983468  0.3683010 0.2406310
## X12 -0.3430122 -0.3647332 0.2824578
# Post
psych::describeBy(socialconnect_scale_include2$EMA_social_comp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed     mad min max
## X11    1      0    1 667 13.09669 4.226561     13 13.25728 4.44780  -2  22
## X12    2      1    1 678 12.25745 4.798022     13 12.51417 4.76008  -2  22
##     range       skew   kurtosis        se
## X11    24 -0.4221188  0.3060899 0.1636531
## X12    24 -0.4717985 -0.1965581 0.1842669
boxplot(EMA_social_comp~Condition,socialconnect_scale_include2, main="Brief Social Connection Scores by Condition",
   xlab="Condition", ylab="Social Connection Comp Score")

Mixed model for Social Connection Composite (POST)

# NEED TO UPDATE THESE TO REMOVE "onewkpost_14wkpost" FROM (1 + onewkpost_14wkpost | PID) once I confirm with David. 

# Model 1 Social Comp - this is a group by time interaction with 3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)

socialconnect_scale_include2$onewkpost_14wkpost <- relevel(socialconnect_scale_include2$onewkpost_14wkpost, ref = "1week")

EQ_MM_Social_model <- lmer (EMA_social_comp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) # "1 + onewkpost_14wkpost | PID" estimates a random intercept and a random slope for each individual


EQ_MM_Social_model
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: 
## EMA_social_comp ~ onewkpost_14wkpost * Condition + (1 + onewkpost_14wkpost |  
##     PID)
##    Data: socialconnect_scale_include2
## REML criterion at convergence: 7280.51
## Random effects:
##  Groups   Name                     Std.Dev. Corr
##  PID      (Intercept)              2.692        
##           onewkpost_14wkpost14week 1.511    0.37
##  Residual                          3.229        
## Number of obs: 1345, groups:  PID, 111
## Fixed Effects:
##                        (Intercept)            onewkpost_14wkpost14week  
##                            13.6231                             -0.9319  
##                          Condition  onewkpost_14wkpost14week:Condition  
##                            -0.6194                             -0.3074
summary(EQ_MM_Social_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## EMA_social_comp ~ onewkpost_14wkpost * Condition + (1 + onewkpost_14wkpost |  
##     PID)
##    Data: socialconnect_scale_include2
## 
## REML criterion at convergence: 7280.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3066 -0.5551  0.0798  0.6238  3.3182 
## 
## Random effects:
##  Groups   Name                     Variance Std.Dev. Corr
##  PID      (Intercept)               7.245   2.692        
##           onewkpost_14wkpost14week  2.284   1.511    0.37
##  Residual                          10.426   3.229        
## Number of obs: 1345, groups:  PID, 111
## 
## Fixed effects:
##                                    Estimate Std. Error       df t value
## (Intercept)                         13.6231     0.4060 110.2270  33.551
## onewkpost_14wkpost14week            -0.9319     0.3251 103.6763  -2.867
## Condition                           -0.6194     0.5704 109.1143  -1.086
## onewkpost_14wkpost14week:Condition  -0.3074     0.4598 103.2911  -0.669
##                                    Pr(>|t|)    
## (Intercept)                         < 2e-16 ***
## onewkpost_14wkpost14week            0.00502 ** 
## Condition                           0.27990    
## onewkpost_14wkpost14week:Condition  0.50522    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) on_1414 Condtn
## onwkps_1414 -0.045               
## Condition   -0.712  0.032        
## onwk_1414:C  0.032 -0.707  -0.038
library(sjPlot)
sjPlot::tab_model(EQ_MM_Social_model)
  EMA social comp
Predictors Estimates CI p
(Intercept) 13.62 12.83 – 14.42 <0.001
onewkpost_14wkpost
[14week] 		-0.93 -1.57 – -0.29 0.004
Condition -0.62 -1.74 – 0.50 0.278
onewkpost_14wkpost
[14week] * Condition 		-0.31 -1.21 – 0.59 0.504
Random Effects
σ2 10.43
τ00 PID 7.24
τ11 PID.onewkpost_14wkpost14week 2.28
ρ01 PID 0.37
ICC 0.49
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.021 / 0.498 
#This model does not yet account for person-level predictors 

# Plotting model

library(sjPlot)
plot_model(model = EQ_MM_Social_model, type = "int",   se = FALSE) #se = FALSE to confirm that we are plotting confidence intervals, not standard error

WRITE UP - Mixed Model Post Social Composite

It was predicted that there would be decreases in social connection over the course of the semester, but that there would be protective effects on social connection decreases in our treatment group relative to our active control group over time. Mixed effect linear models showed that there was no condition main effect (t(109) = -1.09, p = 0.28), but there was a significant main effect of time (t(103) = -2.87, p =0.01). Critically, there was not a significant time by condition interaction (t(103) = -0.67, p =0.51, d = 0.55). Specifically, group 1 had similar decreases from 1-week post (M = 13.6, SE = 0.41) to end of semester (M = 12.7, SE = 0.51) that were comparable to group 2 at 1-week post (M = 13.0, SE = 0.40) to end of semester (M = 11.8, SE = 0.51).

Modeling approach using a dummy coding

Visualize cell means and marginal means

library(apaTables)

soc_connect_means <- apa.2way.table(iv1 = Condition, 
               iv2 = onewkpost_14wkpost, 
               dv = EMA_social_comp, 
               data = socialconnect_scale_include2,
               show.marginal.means = TRUE)

soc_connect_means
## 
## 
## Means and standard deviations for EMA_social_comp as a function of a 2(Condition) X 2(onewkpost_14wkpost) design 
## 
##            onewkpost_14wkpost                               
##                         1week      14week      Marginal     
##  Condition                  M   SD      M   SD        M   SD
##          0              13.57 3.91  12.65 4.46    13.10 4.23
##          1              12.90 4.36  11.58 5.14    12.26 4.80
##   Marginal              13.23 4.16  12.13 4.83              
## 
## Note. M and SD represent mean and standard deviation, respectively. 
## Marginal indicates the means and standard deviations pertaining to main effects.

Creating dummy codes to use in the model (See reference note 1.a).

#Dummy code using if/else statement (for Condition variable): 
socialconnect_scale_include2$Dummy_condition <- if_else(socialconnect_scale_include2$Condition == 1, 1, 0) #Dummy code Condition 1 as "1" and Condition 0 as "0"

#Dummy code using if/else statement (for Time variable):
socialconnect_scale_include2$Dummy_week <- if_else(socialconnect_scale_include2$onewkpost_14wkpost == "14week", 1, 0) #Dummy code 1-week post as "0" and 14-weeks post as "1"

#Checking 
table(socialconnect_scale_include2$Dummy_condition,socialconnect_scale_include2$Dummy_week)
##    
##       0   1
##   0 323 344
##   1 347 331

#Running the multi-level model using dummy-coded variables (see Reference note 1.c).

# Modeling for 1-week post 
model_soc_Dummy_1wkPost <- lmer (EMA_social_comp ~  Dummy_condition*Dummy_week +  (1 | PID), #(1 | PID) accounts for dependency
             data = socialconnect_scale_include2)
summary(model_soc_Dummy_1wkPost) #Gives same output as if Condition were not dummy coded.  
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EMA_social_comp ~ Dummy_condition * Dummy_week + (1 | PID)
##    Data: socialconnect_scale_include2
## 
## REML criterion at convergence: 7307.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4625 -0.5609  0.0320  0.6458  3.1686 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PID      (Intercept)  9.253   3.042   
##  Residual             11.017   3.319   
## Number of obs: 1345, groups:  PID, 111
## 
## Fixed effects:
##                             Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                  13.6226     0.4507  130.1090  30.225  < 2e-16 ***
## Dummy_condition              -0.5956     0.6334  129.0660  -0.940 0.348826    
## Dummy_week                   -0.9290     0.2590 1235.7642  -3.587 0.000348 ***
## Dummy_condition:Dummy_week   -0.3506     0.3659 1237.9493  -0.958 0.338120    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dmmy_c Dmmy_w
## Dummy_cndtn -0.712              
## Dummy_week  -0.297  0.211       
## Dmmy_cnd:D_  0.210 -0.289 -0.708
library(sjPlot)
sjPlot::tab_model(model_soc_Dummy_1wkPost)
  EMA social comp
Predictors Estimates CI p
(Intercept) 13.62 12.74 – 14.51 <0.001
Dummy_condition -0.60 -1.84 – 0.65 0.347
Dummy_week -0.93 -1.44 – -0.42 <0.001
Dummy_condition *
Dummy_week 		-0.35 -1.07 – 0.37 0.338
Random Effects
σ2 11.02
τ00 PID 9.25
ICC 0.46
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.022 / 0.468 
plot_model(model_soc_Dummy_1wkPost, type = "pred", terms = c("Dummy_week[1,0]", "Dummy_condition"))

Using dummy codes to test the null hypotheses that there is no significant difference between C0 and C1 for social connection at 1-week post

#fit <- lm(model_soc_Dummy_1wkPost) #Looking at the lm output (instead of anova(model)) for it (because ANOVAs are just fancy linear models) 
#summary(fit)
#lm(model_soc_Dummy_1wkPost) %>% 
  #broom::tidy()

#Interpreting the output B0 (intercept) + B1 (Dummy_condition) * onewkpost_14wkpost (Dummy_week) Intercept (B0): The mean for Condition 0 at 1-week post (i.e. when both dummys are zero)
Dummy_week (B1): The difference across time for Condition 0 Dummy_condition (simple effect - the main effect of condition): Codes for the difference between the groups at 1-week post. Here, p = 0.06585

Using dummy codes to test the null hypotheses that there is no significant difference between C0 and C1 for social connection at 14-weeks post

#OPTION

#Dummy code using if/else statement (for Time variable):
socialconnect_scale_include2$Dummy_week <- if_else(socialconnect_scale_include2$onewkpost_14wkpost == "1week", 1, 0) #Dummy code 1-week post as "1" and 14-weeks post as "0"

#Checking 
#table(socialconnect_scale_include2$Dummy_condition,socialconnect_scale_include2$Dummy_week)

# Modeling for 14-week post 
model_soc_Dummy_14wkPost <- lmer (EMA_social_comp ~  Dummy_condition*Dummy_week +  (1 | PID), #(1 | PID) accounts for dependency
             data = socialconnect_scale_include2)
summary(model_soc_Dummy_14wkPost) #Gives same output as if Condition were not dummy coded.  
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EMA_social_comp ~ Dummy_condition * Dummy_week + (1 | PID)
##    Data: socialconnect_scale_include2
## 
## REML criterion at convergence: 7307.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4625 -0.5609  0.0320  0.6458  3.1686 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PID      (Intercept)  9.253   3.042   
##  Residual             11.017   3.319   
## Number of obs: 1345, groups:  PID, 111
## 
## Fixed effects:
##                             Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                  12.6936     0.4482  127.2565  28.324  < 2e-16 ***
## Dummy_condition              -0.9461     0.6333  128.9657  -1.494 0.137607    
## Dummy_week                    0.9290     0.2590 1235.7642   3.587 0.000348 ***
## Dummy_condition:Dummy_week    0.3506     0.3659 1237.9493   0.958 0.338120    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dmmy_c Dmmy_w
## Dummy_cndtn -0.708              
## Dummy_week  -0.279  0.198       
## Dmmy_cnd:D_  0.198 -0.289 -0.708
library(sjPlot)
sjPlot::tab_model(model_soc_Dummy_14wkPost)
  EMA social comp
Predictors Estimates CI p
(Intercept) 12.69 11.82 – 13.57 <0.001
Dummy_condition -0.95 -2.19 – 0.30 0.135
Dummy_week 0.93 0.42 – 1.44 <0.001
Dummy_condition *
Dummy_week 		0.35 -0.37 – 1.07 0.338
Random Effects
σ2 11.02
τ00 PID 9.25
ICC 0.46
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.022 / 0.468 
#fit_14 <- lm(model_soc_Dummy_14wkPost) #Looking at the lm output (instead of anova(model)) for it (because ANOVAs are just fancy linear models) 
#summary(fit_14)
#lm(model_soc_Dummy_14wkPost) %>% 
  #broom::tidy()

Reference Notes For Dummy Coding

#1.a Dummy codes turn things into zeros and ones - Condition is already in 0s and 1s, BUT now the Condition variable codes for just one of the groups (with respect to the reference group). The significance test of a dummy code captures the difference between the coded group and the refernce group. So here, with two groups, it will code for the difference between the two groups. For this 2x2 (group x time), the main effect of the dummy code will actually be coding for the main effect of that variable (e.g. group) at the zero value of the other variable of the other dummy code (e.g. time) - so essentially it tests the effect of Condition at time 1.

#1.b Creating dummy code using dummy.code from psych package (https://www.rdocumentation.org/packages/psych/versions/2.0.8/topics/dummy.code)
  # this is not nested within dplyr. 
  # Condition is already zeros and ones, so they will just be put in the regression model.

#1.c Resource: Stats lecture Interactions (III) slide 24 and Interactions (IV) slides 11-14. This is basically a factorial ANOVA - i.e. 2x2 matrix of variable 1 (timepoint -oneweek versus fourteen weeks) and variable 2 (Condition: 0 versus 1), but it is a multi-level model because people are in more than 1 of those boxes because they have multiple timepoints. This is a conditional growth model - dropping random slope (e.g. https://rpsychologist.com/r-guide-longitudinal-lme-lmer)
  # Note: Estimates and SEs expected to be similar (but not identical to) breaking the data apart and running a simple regression.

EXPLORING

# Marginal Means

library(emmeans)
## Warning: package 'emmeans' was built under R version 3.6.2
soc_rg <- ref_grid(model_soc_Dummy_1wkPost)
summary(soc_rg) #gives the  and the mean + SE of C0 and C1 at 14-week post
##  Dummy_condition Dummy_week prediction    SE  df
##                0          0       13.6 0.451 131
##                1          0       13.0 0.445 128
##                0          1       12.7 0.448 128
##                1          1       11.7 0.447 131
## 
## Degrees-of-freedom method: kenward-roger
# Obtaining marginal means (along with confidence intervals) to interpret main effects. Note: Interpretation of the main effects and interaction in a factorial design will usually require follow-up comparisons conducted at the level of the effect.

timepoint_lsm = emmeans::lsmeans(soc_rg, "Dummy_condition")
## NOTE: Results may be misleading due to involvement in interactions
timepoint_lsm
##  Dummy_condition lsmean    SE  df lower.CL upper.CL
##                0   13.2 0.430 109     12.3     14.0
##                1   12.4 0.427 109     11.5     13.2
## 
## Results are averaged over the levels of: Dummy_week 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
pairs(timepoint_lsm, adjust = "holm") #Significant Main Effect of Time (for C1?) 
##  contrast estimate    SE  df t.ratio p.value
##  0 - 1       0.771 0.606 109 1.271   0.2063 
## 
## Results are averaged over the levels of: Dummy_week 
## Degrees-of-freedom method: kenward-roger

Effect Size

#Note From Elliot - effect sizes in an MLM do not need to be cohen's d. You can report the standardized regression coefficients (the beta) in standardized units. e.g these would tell you the simple effect of time pre to post - i.e. the change in one of the groups (or both of the groups together). OR, you could get (at 1-week post and at 14-weeks post) what was the difference between the groups. Get this in terms of a standardized b. 

# Reporting effect size in standardized regression coefficients (i.e. betas) using std_beta function from here: https://strengejacke.github.io/sjstats/reference/std_beta.html, which can take fitted linear (mixed) model of class lm, merMod (lme4 package), gls or stanreg. Note: Standardized coefficients refer to how many standard deviations a dependent variable will change, per standard deviation increase in the predictor variable. Standardization of the coefficient is usually done to answer the question of which of the independent variables have a greater effect on the dependent variable in a multiple regression analysis, when the variables are measured in different units of measurement (for example, income measured in dollars and family size measured in number of individuals). Standardized regression coefficients (i.e. betas) compare the relative importance of each coefficient (e.g. time or condition) in a regression model.
#Interpretation of effect sizes, as a rule-of-thumb:
  # small effect >= 0.1
  # medium effect >= 0.3
  # large effect >= 0.5

# beta for time, group, and interaction
#install.packages("sjstats")
sjstats::std_beta(EQ_MM_Social_model, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients. 
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.14031010
## 2           onewkpost_14wkpost14week     -0.23931613
## 3                          Condition     -0.06821496
## 4 onewkpost_14wkpost14week:Condition     -0.03385820
#Below code is OLD - (not using r-squared)
# Attempt #1

  # Calculating Cohen's d (effect size) - using method from Jarrett Byrnes described here:     https://stats.stackexchange.com/questions/95054/how-to-get-an-overall-p-value-and-effect-size-for-a-categorical-factor-in-a-mi

#r2.corr.mer <- function(m) {
  #lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))
  #summary(lmfit)$r.squared
#}

#r2.corr.mer(EQ_MM_Social_model)

  # Attempt #2
#install.packages("MuMIn")

#MuMIn::r.squaredGLMM(EQ_MM_Social_model)  
# 'R2m' = marginal R squared (proportion of variance explained by the fixed factor(s) alone)
# 'R2c' = conditional R squared (proportion of variance explained by both the fixed and random factors)

#Unsuccessful attempts for reference
#lsr::cohensD(x      = socialconnect_scale_include2$EMA_social_comp_Condition0,
            # y      = socialconnect_scale_include2$EMA_social_comp_Condition1, 
             # method = "paired")

#Option from Stats lab
#lsr::etaSquared(model_soc_dummy)

Reporting the Effect size using the standardized regression coefficients (the beta) - outliers included

  1. From 1-week post to 14-week post (aka an increase of one standard deviation in time), social connection decreases by 0.23 standard deviations. (“MAIN EFFECT OF TIME”)
  2. From condition 1 to condition 2 (aka an increase of one standard deviation in condition (1 or 2)), social connection decreases by 0.07 standard deviations.
  3. The overall effect size (standardized regression coefficient) for the interaction between time and condition was -0.03 (in standardized units). #QUESTION: How do I intertpret this? Write this up?

Plotting EMA_comp scores

library(ggplot2)

#+++++++++++++++++++++++++
# Function to calculate the mean and the standard deviation for each group
  ## Summarizes data - Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) { 
  ##   data: a data frame.
  ##   measurevar: the name of a column that contains the variable to be summariezed
  ##   groupvars: a vector containing names of columns that contain grouping variables
  ##   na.rm: a boolean that indicates whether to ignore NA's
  ##   conf.interval: the percent range of the confidence interval (default is 95%)
    library(plyr)
    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }
   ## This does the summary. For each group's data frame, return a vector with N, mean, and sd
    datac <- ddply(data, groupvars, .drop=.drop,
      .fun = function(xx, col) {
        c(N    = length2(xx[[col]], na.rm=na.rm),
          mean = mean   (xx[[col]], na.rm=na.rm),
          sd   = sd     (xx[[col]], na.rm=na.rm)
        )
      },
      measurevar
    )
    # Rename the "mean" column    
    datac <- rename(datac, c("mean" = measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}

#+++++++++++++++++++++++++

soc_connect_summary <- summarySE(socialconnect_scale_include2, measurevar="EMA_social_comp", groupvars=c("onewkpost_14wkpost", "Condition"))
soc_connect_summary
##   onewkpost_14wkpost Condition   N EMA_social_comp       sd        se        ci
## 1              1week         0 323        13.57223 3.910311 0.2175755 0.4280490
## 2              1week         1 347        12.90287 4.359070 0.2340071 0.4602555
## 3             14week         0 344        12.65017 4.463040 0.2406310 0.4732981
## 4             14week         1 331        11.58084 5.138869 0.2824578 0.5556450
soc_connect_summary$Condition <- as.factor(soc_connect_summary$Condition)

library(plyr)
soc_connect_summary$onewkpost_14wkpost<- revalue(soc_connect_summary$onewkpost_14wkpost, c("1week"="1-Week Post", "14week"="14-Weeks Post"))

#Basic line plot with confidence intervals
p<- ggplot(data=soc_connect_summary, aes(x=as.factor(onewkpost_14wkpost), y=EMA_social_comp, group=Condition, color=Condition)) +
  geom_line(aes(linetype=Condition)) + geom_point(aes(shape=Condition))+
  scale_color_brewer(palette="Paired")+
  geom_point(position=position_dodge(0.05), size=3, shape=21, fill="white") + # 21 is filled circle
  labs(title="Social Connection at 1-week and 14-weeks Post",x="Timepoint", y = "Social Connection")+
  theme_classic() +
  geom_errorbar(aes(ymin=EMA_social_comp-ci, ymax=EMA_social_comp+ci), width=.2,
                 position=position_dodge(0.05)) # The errorbars overlapped, so use position_dodge to move them horizontally .05 to the left and right
p + theme_bw() + scale_color_manual(values=c('orange1','olivedrab3')) +
    theme(legend.justification=c(1,0),
          legend.position=c(1,0)) 

# QUESTION: Does it make sense to have an N? Do I use N for Pps or for observations
#QUESTION: Why are the error bars (which should be confidence intervals, not overlapping at time 2 - as in the sjplot model graph?)

Removing Outliers

#Obtaining values of outliers 

library(dplyr)
library(ggplot2)

is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

socialconnect_scale_include2 %>%
  group_by(Condition) %>%
  mutate(outlier = ifelse(is_outlier(EMA_social_comp), EMA_social_comp, as.numeric(NA))) %>%
  ggplot(., aes(x = factor(Condition), y = EMA_social_comp)) +
    geom_boxplot() +
    geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)

# NEED TO FIGURE OUT HOW TO Remove outliers at both levels

Removing Outliers

socialconnect_scale_include2_out <- socialconnect_scale_include2 %>% 
  filter(Condition == 0, EMA_social_comp > 3) 

socialconnect_scale_include2_1 <- socialconnect_scale_include2 %>% 
  filter(Condition == 1, EMA_social_comp > -2) 

socialconnect_scale_out <- rbind(socialconnect_scale_include2_out, socialconnect_scale_include2_1)
#run the models with and without these outliers

Mixed model for EMA Social Connection (outliers removed)

library(lme4)
library(lmerTest)

socialconnect_scale_out <- socialconnect_scale_out %>% 
  mutate(onewkpost_14wkpost = as.factor(onewkpost_14wkpost), 
         timepoint = as.factor(timepoint))

EQ_MM_SocialConnect_model_out <- lmer (EMA_social_comp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_out) 
 
#this is a group by time interaction
#3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)
# "1 + onewkpost_14wkpost | PID" estimate a random intercept and a random slope for each individual 
summary(EQ_MM_SocialConnect_model_out)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## EMA_social_comp ~ onewkpost_14wkpost * Condition + (1 + onewkpost_14wkpost |  
##     PID)
##    Data: socialconnect_scale_out
## 
## REML criterion at convergence: 7098.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3786 -0.5609  0.0617  0.6388  3.4212 
## 
## Random effects:
##  Groups   Name                     Variance Std.Dev. Corr
##  PID      (Intercept)              7.058    2.657        
##           onewkpost_14wkpost14week 1.919    1.385    0.29
##  Residual                          9.736    3.120        
## Number of obs: 1328, groups:  PID, 111
## 
## Fixed effects:
##                                    Estimate Std. Error       df t value
## (Intercept)                         13.6835     0.3994 110.0831  34.263
## onewkpost_14wkpost14week            -0.7267     0.3102 104.0786  -2.343
## Condition                           -0.6728     0.5609 108.8821  -1.200
## onewkpost_14wkpost14week:Condition  -0.4486     0.4374 102.3356  -1.026
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## onewkpost_14wkpost14week              0.021 *  
## Condition                             0.233    
## onewkpost_14wkpost14week:Condition    0.307    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) on_1414 Condtn
## onwkps_1414 -0.097               
## Condition   -0.712  0.069        
## onwk_1414:C  0.069 -0.709  -0.089
## convergence code: 0
## Model failed to converge with max|grad| = 0.00379648 (tol = 0.002, component 1)
sjPlot::tab_model(EQ_MM_SocialConnect_model_out)
  EMA social comp
Predictors Estimates CI p
(Intercept) 13.68 12.90 – 14.47 <0.001
onewkpost_14wkpost
[14week] 		-0.73 -1.33 – -0.12 0.019
Condition -0.67 -1.77 – 0.43 0.230
onewkpost_14wkpost
[14week] * Condition 		-0.45 -1.31 – 0.41 0.305
Random Effects
σ2 9.74
τ00 PID 7.06
τ11 PID.onewkpost_14wkpost14week 1.92
ρ01 PID 0.29
ICC 0.48
N PID 111
Observations 1328
Marginal R2 / Conditional R2 0.022 / 0.494 
#This model does not yet account for person-level predictors

Calculating the effect size for group, time, and interaction from the model

# Calculating Cohen's d (effect size) - using method from Jarrett Byrnes described here:     https://stats.stackexchange.com/questions/95054/how-to-get-an-overall-p-value-and-effect-size-for-a-categorical-factor-in-a-mi

# Attempt #1 - not using r-squared (see note below)
#r2.corr.mer(EQ_MM_SocialConnect_model_out) #get the R-squared (i.e. the total amount of variance explained by a model), which is an effect size (Note: but the variance explained by the model could be the time, the condition, or the interaction). 

# Attempt #2 - not using r-squared (see note below)
#MuMIn::r.squaredGLMM(EQ_MM_SocialConnect_model_out)  
# 'R2m' = marginal R squared (proportion of variance explained by the fixed factor(s) alone)
# 'R2c' = conditional R squared (proportion of variance explained by both the fixed and random factors)

#Note From Elliot - effect sizes in an MLM do not need to be cohen's d. You can report the standardized regression coefficients (the beta) in standardized units. e.g these would tell you the simple effect of time pre to post - i.e. the change in one of the groups (or both of the groups together). OR, you could get (at 1-week post and at 14-weeks post) what was the difference between the groups. Get this in terms of a standardized b. 

# Reporting effect size in standardized regression coefficients (i.e. betas) using std_beta function from here: https://strengejacke.github.io/sjstats/reference/std_beta.html, which can take fitted linear (mixed) model of class lm, merMod (lme4 package), gls or stanreg. Note: Standardized coefficients refer to how many standard deviations a dependent variable will change, per standard deviation increase in the predictor variable. Standardization of the coefficient is usually done to answer the question of which of the independent variables have a greater effect on the dependent variable in a multiple regression analysis, when the variables are measured in different units of measurement (for example, income measured in dollars and family size measured in number of individuals). Standardized regression coefficients (i.e. betas) compare the relative importance of each coefficient (e.g. time or condition) in a regression model. 

#Find a package for a function
#install.packages("sos")
#require("sos")
#findFn("std_beta")

# beta for time, group, and interaction
#install.packages("sjstats")
sjstats::std_beta(EQ_MM_SocialConnect_model_out, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients. 
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.11802018
## 2           onewkpost_14wkpost14week     -0.21926308
## 3                          Condition     -0.07726937
## 4 onewkpost_14wkpost14week:Condition     -0.05152899

Reporting the Effect size using the standardized regression coefficients (the beta) - outliers removed

With every increase of one standard deviation in time (1-week post to 14-week post), social connection decreases by 0.22 standard deviations. With every increase of one standard deviation in condition (1 or 2), social connection decreases by 0.08 standard deviations. The overall effect size (standardized regression coefficient) for the interaction between time and condition was -0.05 (in standardized units). #QUESTION: How do I intertpret this? Write this up?

Modeling approach using a dummy coding (outliers removed)

Visualize cell means and marginal means

library(apaTables)

soc_connect_means_out <- apa.2way.table(iv1 = Condition, 
               iv2 = onewkpost_14wkpost, 
               dv = EMA_social_comp, 
               data = socialconnect_scale_out,
               show.marginal.means = TRUE)

soc_connect_means_out
## 
## 
## Means and standard deviations for EMA_social_comp as a function of a 2(Condition) X 2(onewkpost_14wkpost) design 
## 
##            onewkpost_14wkpost                               
##                         1week      14week      Marginal     
##  Condition                  M   SD      M   SD        M   SD
##          0              13.65 3.78  13.06 3.96    13.35 3.88
##          1              12.90 4.36  11.71 4.99    12.32 4.71
##   Marginal              13.26 4.11  12.39 4.55              
## 
## Note. M and SD represent mean and standard deviation, respectively. 
## Marginal indicates the means and standard deviations pertaining to main effects.

Creating dummy codes to use in the model (See reference note 1.a).

#Dummy code using if/else statement (for Condition variable): 
socialconnect_scale_out$Dummy_condition <- if_else(socialconnect_scale_out$Condition == 1, 1, 0) #Dummy code Condition 1 as "1" and Condition 0 as "0"

#Dummy code using if/else statement (for Time variable):
socialconnect_scale_out$Dummy_week <- if_else(socialconnect_scale_out$onewkpost_14wkpost == "14week", 1, 0) #Dummy code 1-week post as "0" and 14-weeks post as "1"

#Checking 
table(socialconnect_scale_out$Dummy_condition,socialconnect_scale_out$Dummy_week)
##    
##       0   1
##   0 321 332
##   1 347 328

#Running the multi-level model using dummy-coded variables (see Reference note 1.c).

# Modeling for 1-week post 
model_soc_Dummy_1wkPost_out <- lmer (EMA_social_comp ~  Dummy_condition*Dummy_week +  (1 | PID), #(1 | PID) accounts for dependency
             data = socialconnect_scale_out)
summary(model_soc_Dummy_1wkPost_out) #Gives same output as if Condition were not dummy coded.  
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EMA_social_comp ~ Dummy_condition * Dummy_week + (1 | PID)
##    Data: socialconnect_scale_out
## 
## REML criterion at convergence: 7118.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4938 -0.5921  0.0256  0.6608  3.2754 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PID      (Intercept)  8.493   2.914   
##  Residual             10.234   3.199   
## Number of obs: 1328, groups:  PID, 111
## 
## Fixed effects:
##                             Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                  13.6801     0.4325  129.8190  31.632  < 2e-16 ***
## Dummy_condition              -0.6468     0.6076  128.6307  -1.065  0.28905    
## Dummy_week                   -0.7051     0.2526 1219.2220  -2.791  0.00533 ** 
## Dummy_condition:Dummy_week   -0.5058     0.3551 1220.8975  -1.424  0.15463    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dmmy_c Dmmy_w
## Dummy_cndtn -0.712              
## Dummy_week  -0.297  0.212       
## Dmmy_cnd:D_  0.212 -0.290 -0.711
library(sjPlot)
sjPlot::tab_model(model_soc_Dummy_1wkPost_out)
  EMA social comp
Predictors Estimates CI p
(Intercept) 13.68 12.83 – 14.53 <0.001
Dummy_condition -0.65 -1.84 – 0.54 0.287
Dummy_week -0.71 -1.20 – -0.21 0.005
Dummy_condition *
Dummy_week 		-0.51 -1.20 – 0.19 0.154
Random Effects
σ2 10.23
τ00 PID 8.49
ICC 0.45
N PID 111
Observations 1328
Marginal R2 / Conditional R2 0.023 / 0.466 
plot_model(model_soc_Dummy_1wkPost_out, type = "pred", terms = c("Dummy_week[1,0]", "Dummy_condition"))

Using dummy codes to test the null hypotheses that there is no significant difference between C0 and C1 for social connection at 1-week post

fit_out <- lm(model_soc_Dummy_1wkPost_out) #Looking at the lm output (instead of anova(model)) for it (because ANOVAs are just fancy linear models) 
summary(fit_out)
## 
## Call:
## lm(formula = model_soc_Dummy_1wkPost_out)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.8002  -2.7391   0.1976   3.2044  10.5986 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                13.128609   0.402848  32.589   <2e-16 ***
## Dummy_condition            -0.726639   0.333316  -2.180   0.0294 *  
## Dummy_week                 -0.589888   0.336590  -1.753   0.0799 .  
## PID                         0.001058   0.000652   1.623   0.1049    
## Dummy_condition:Dummy_week -0.609619   0.472174  -1.291   0.1969    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.3 on 1323 degrees of freedom
## Multiple R-squared:  0.02785,    Adjusted R-squared:  0.02491 
## F-statistic: 9.477 on 4 and 1323 DF,  p-value: 1.489e-07
lm(model_soc_Dummy_1wkPost_out) %>% 
  broom::tidy()
## # A tibble: 5 x 5
##   term                       estimate std.error statistic   p.value
##   <chr>                         <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)                13.1      0.403        32.6  1.62e-171
## 2 Dummy_condition            -0.727    0.333        -2.18 2.94e-  2
## 3 Dummy_week                 -0.590    0.337        -1.75 7.99e-  2
## 4 PID                         0.00106  0.000652      1.62 1.05e-  1
## 5 Dummy_condition:Dummy_week -0.610    0.472        -1.29 1.97e-  1

#Interpreting the output B0 (intercept) + B1 (Dummy_condition) * onewkpost_14wkpost (Dummy_week) Intercept (B0): The mean for Condition 0 at 1-week post (i.e. when both dummys are zero)
Dummy_week (B1): The difference across time for Condition 0 Dummy_condition (simple effect - the main effect of condition): Codes for the difference between the groups at 1-week post.

Using dummy codes to test the null hypotheses that there is no significant difference between C0 and C1 for social connection at 14-weeks post

#Dummy code using if/else statement (for Time variable):
socialconnect_scale_out$Dummy_week <- if_else(socialconnect_scale_out$onewkpost_14wkpost == "1week", 1, 0) #Dummy code 1-week post as "1" and 14-weeks post as "0"

#Checking 
#table(socialconnect_scale_out$Dummy_condition,socialconnect_scale_out$Dummy_week)

# Modeling for 14-week post 
model_soc_Dummy_14wkPost_out <- lmer (EMA_social_comp ~  Dummy_condition*Dummy_week +  (1 | PID), #(1 | PID) accounts for dependency
             data = socialconnect_scale_out)
summary(model_soc_Dummy_14wkPost_out) #Gives same output as if Condition were not dummy coded.  
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EMA_social_comp ~ Dummy_condition * Dummy_week + (1 | PID)
##    Data: socialconnect_scale_out
## 
## REML criterion at convergence: 7118.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4938 -0.5921  0.0256  0.6608  3.2754 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PID      (Intercept)  8.493   2.914   
##  Residual             10.234   3.199   
## Number of obs: 1328, groups:  PID, 111
## 
## Fixed effects:
##                             Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                  12.9750     0.4311  128.2423  30.096  < 2e-16 ***
## Dummy_condition              -1.1526     0.6085  129.4020  -1.894  0.06045 .  
## Dummy_week                    0.7051     0.2526 1219.2220   2.791  0.00533 ** 
## Dummy_condition:Dummy_week    0.5058     0.3551 1220.8975   1.424  0.15463    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dmmy_c Dmmy_w
## Dummy_cndtn -0.708              
## Dummy_week  -0.288  0.204       
## Dmmy_cnd:D_  0.205 -0.294 -0.711
library(sjPlot)
sjPlot::tab_model(model_soc_Dummy_14wkPost_out)
  EMA social comp
Predictors Estimates CI p
(Intercept) 12.97 12.13 – 13.82 <0.001
Dummy_condition -1.15 -2.35 – 0.04 0.058
Dummy_week 0.71 0.21 – 1.20 0.005
Dummy_condition *
Dummy_week 		0.51 -0.19 – 1.20 0.154
Random Effects
σ2 10.23
τ00 PID 8.49
ICC 0.45
N PID 111
Observations 1328
Marginal R2 / Conditional R2 0.023 / 0.466 
fit_14_out <- lm(model_soc_Dummy_14wkPost_out) #Looking at the lm output (instead of anova(model)) for it (because ANOVAs are just fancy linear models) 
summary(fit_14_out)
## 
## Call:
## lm(formula = model_soc_Dummy_14wkPost_out)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.8002  -2.7391   0.1976   3.2044  10.5986 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                12.538722   0.398196  31.489  < 2e-16 ***
## Dummy_condition            -1.336258   0.334929  -3.990 6.98e-05 ***
## Dummy_week                  0.589888   0.336590   1.753   0.0799 .  
## PID                         0.001058   0.000652   1.623   0.1049    
## Dummy_condition:Dummy_week  0.609619   0.472174   1.291   0.1969    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.3 on 1323 degrees of freedom
## Multiple R-squared:  0.02785,    Adjusted R-squared:  0.02491 
## F-statistic: 9.477 on 4 and 1323 DF,  p-value: 1.489e-07
lm(model_soc_Dummy_14wkPost_out) %>% 
  broom::tidy()
## # A tibble: 5 x 5
##   term                       estimate std.error statistic   p.value
##   <chr>                         <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)                12.5      0.398        31.5  6.96e-163
## 2 Dummy_condition            -1.34     0.335        -3.99 6.98e-  5
## 3 Dummy_week                  0.590    0.337         1.75 7.99e-  2
## 4 PID                         0.00106  0.000652      1.62 1.05e-  1
## 5 Dummy_condition:Dummy_week  0.610    0.472         1.29 1.97e-  1

Plotting EMA_comp scores w/ outliers removed

library(ggplot2)

soc_connect_summary_out <- summarySE(socialconnect_scale_out, measurevar="EMA_social_comp", groupvars=c("onewkpost_14wkpost", "Condition"))
soc_connect_summary_out
##   onewkpost_14wkpost Condition   N EMA_social_comp       sd        se        ci
## 1              1week         0 321        13.65368 3.782837 0.2111374 0.4153927
## 2              1week         1 347        12.90287 4.359070 0.2340071 0.4602555
## 3             14week         0 332        13.05921 3.963221 0.2175100 0.4278763
## 4             14week         1 328        11.70506 4.994264 0.2757621 0.5424916
soc_connect_summary_out$Condition <- as.factor(soc_connect_summary_out$Condition)

library(plyr)
soc_connect_summary_out$onewkpost_14wkpost<- revalue(soc_connect_summary_out$onewkpost_14wkpost, c("1week"="1-Week Post", "14week"="14-Weeks Post"))

#Basic line plot 
#p<- ggplot(data=soc_connect_summary_out, aes(x=as.factor(onewkpost_14wkpost), y=EMA_social_comp, group=Condition, color=Condition)) +
  #geom_line(aes(linetype=Condition)) + geom_point(aes(shape=Condition))+
  #scale_color_brewer(palette="Paired")+
  #labs(title="Social Connection at 1-week and 14-weeks Post",x="Timepoint", y = "Social Connection")+
  #theme_classic() 
#p + theme_classic() + scale_color_manual(values=c('orange1','olivedrab3'))

#Basic line plot with confidence intervals
p<- ggplot(data=soc_connect_summary_out, aes(x=as.factor(onewkpost_14wkpost), y=EMA_social_comp, group=Condition, color=Condition)) +
  geom_line(aes(linetype=Condition)) + geom_point(aes(shape=Condition))+
  scale_color_brewer(palette="Paired")+
  labs(title="Social Connection at 1-week and 14-weeks Post",x="Timepoint", y = "Social Connection")+
  theme_classic() +
  geom_errorbar(aes(ymin=EMA_social_comp-ci, ymax=EMA_social_comp+ci), width=.2,
                 position=position_dodge(0.05))
p + theme_classic() + scale_color_manual(values=c('orange1','olivedrab3'))

# Does it make sense to have an N? Do I use N for Pps or for observations
#library(sjPlot)
#plot_model(model = EQ_MM_SocialConnect_model_out, type = "int")

Creating separate dataframes for 1-week and 14-week follow-up (outliers removed)

socialconnect_scale_out_1wkpost <- socialconnect_scale_out %>% 
  filter(onewkpost_14wkpost == '1week')

socialconnect_scale_out_14wkpost <- socialconnect_scale_out %>% 
  filter(onewkpost_14wkpost == '14week')

Pulling contrasts (means / SEs) from the model for Social Composite (outliers removed)

Note: emmeans is a post-hoc logical comparison - i.e. you usually do it to decompose a significant interaction. The justification for doing this here is because of the limitations of the data (i.e. not pre measure possible). This is somewhat like a two-sample t-test for each timepoint, but takes the nested structure of the data into account.

///NEGATIVE SOCIAL INTERACTIONS COMPOSITE///

Describing the data

# Descriptives - EMA Social Connection Composite (EMA_social_comp)

# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$EMA_neg_social_comp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd   median  trimmed       mad min
## X11    1      0    1 323 1.956048 1.221739 1.259878 1.740402 0.3852957   1
## X12    2      1    1 347 1.926877 1.099753 2.000000 1.757586 1.4826000   1
##          max    range      skew    kurtosis         se
## X11 6.000000 5.000000 1.2762762  0.94670774 0.06797938
## X12 5.259878 4.259878 0.9822507 -0.01134774 0.05903781
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$EMA_neg_social_comp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed mad min max range
## X11    1      0    1 344 1.871247 1.180213      1 1.645613   0   1   6     5
## X12    2      1    1 331 1.693735 1.132470      1 1.463573   0   1   8     7
##         skew  kurtosis         se
## X11 1.316492 0.9790247 0.06363285
## X12 1.975102 4.5200139 0.06224619
# Post
psych::describeBy(socialconnect_scale_include2$EMA_neg_social_comp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed mad min max range
## X11    1      0    1 667 1.912313 1.200348      1 1.691502   0   1   6     5
## X12    2      1    1 678 1.813057 1.121102      1 1.605171   0   1   8     7
##         skew  kurtosis         se
## X11 1.300701 0.9831728 0.04647765
## X12 1.459353 2.0726548 0.04305567
boxplot(EMA_neg_social_comp~Condition,socialconnect_scale_include2, main="Brief Negative Social Connection Scores by Condition",
   xlab="Condition", ylab="Negative Social Interactions Comp Score")

Mixed model for Negative Social Interactions Composite (POST)

# NEED TO UPDATE THESE TO REMOVE "onewkpost_14wkpost" FROM (1 + onewkpost_14wkpost | PID) once I confirm with David. 

# Model 1 Comp

EQ_MM_Social_model_negcomp <- lmer (EMA_neg_social_comp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) 
 
#this is a group by time interaction
#3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)
# "1 + onewkpost_14wkpost | PID" estimate a random intercept and a random slope for each individual 
#summary(EQ_MM_Social_model)

library(sjPlot)
sjPlot::tab_model(EQ_MM_Social_model_negcomp)
  EMA neg social comp
Predictors Estimates CI p
(Intercept) 1.95 1.74 – 2.16 <0.001
onewkpost_14wkpost
[14week] 		-0.07 -0.28 – 0.13 0.484
Condition 0.00 -0.29 – 0.30 0.993
onewkpost_14wkpost
[14week] * Condition 		-0.15 -0.44 – 0.14 0.310
Random Effects
σ2 0.89
τ00 PID 0.47
τ11 PID.onewkpost_14wkpost14week 0.30
ρ01 PID -0.42
ICC 0.34
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.006 / 0.348 
#This model does not yet account for person-level predictors 

# Plotting model

#library(sjPlot)
#plot_model(model = EQ_MM_Social_model_negcomp, type = "int")

#Using better plot below

# Calculating Cohen's d (effect size) - using method from Jarrett Byrnes described here:     https://stats.stackexchange.com/questions/95054/how-to-get-an-overall-p-value-and-effect-size-for-a-categorical-factor-in-a-mi

# Attempt #1
#r2.corr.mer(EQ_MM_Social_model_negcomp)

# Attempt #2
#MuMIn::r.squaredGLMM(EQ_MM_Social_model_negcomp)  
# 'R2m' = marginal R squared (proportion of variance explained by the fixed factor(s) alone)
# 'R2c' = conditional R squared (proportion of variance explained by both the fixed and random factors)

# Standardized effect size 
sjstats::std_beta(EQ_MM_Social_model_negcomp, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients. 
##                            Parameter Std_Coefficient
## 1                        (Intercept)    0.0740218419
## 2           onewkpost_14wkpost14week   -0.1288100614
## 3                          Condition    0.0006049417
## 4 onewkpost_14wkpost14week:Condition   -0.0649484143

It was predicted that there would be decreases in social connection over the course of the semester, but that there would be protective effects on social connection decreases in our treatment group relative to our active control group over time. Mixed effect linear models showed that there was no condition main effect (t(109) = -1.09, p = 0.28), but there was a significant main effect of time (t(103) = -2.87, p =0.01). Critically, there was not a significant time by condition interaction (t(103) = -0.67, p =0.51, d = 0.55). Specifically, group 1 had similar decreases from 1-week post (M = 13.6, SE = 0.41) to end of semester (M = 12.7, SE = 0.51) that were comparable to group 2 at 1-week post (M = 13.0, SE = 0.40) to end of semester (M = 11.8, SE = 0.51).

Plotting EMA_neg_social_comp scores

library(ggplot2)

neg_soc_connect_summary <- summarySE(socialconnect_scale_include2, measurevar="EMA_neg_social_comp", groupvars=c("onewkpost_14wkpost", "Condition"))
neg_soc_connect_summary
##   onewkpost_14wkpost Condition   N EMA_neg_social_comp       sd         se
## 1              1week         0 323            1.956048 1.221739 0.06797938
## 2              1week         1 347            1.926877 1.099753 0.05903781
## 3             14week         0 344            1.871247 1.180213 0.06363285
## 4             14week         1 331            1.693735 1.132470 0.06224619
##          ci
## 1 0.1337398
## 2 0.1161182
## 3 0.1251597
## 4 0.1224494
neg_soc_connect_summary$Condition <- as.factor(neg_soc_connect_summary$Condition)

library(plyr)
neg_soc_connect_summary$onewkpost_14wkpost<- revalue(neg_soc_connect_summary$onewkpost_14wkpost, c("1week"="1-Week Post", "14week"="14-Weeks Post"))

#Basic line plot with confidence intervals
p<- ggplot(data=neg_soc_connect_summary, aes(x=as.factor(onewkpost_14wkpost), y=EMA_neg_social_comp, group=Condition, color=Condition)) +
  geom_line(aes(linetype=Condition)) + geom_point(aes(shape=Condition))+
  scale_color_brewer(palette="Paired")+
  labs(title="Negative Social Interactions at 1-week and 14-weeks Post",x="Timepoint", y = "Negative Social Interactions")+
  theme_classic() +
  geom_errorbar(aes(ymin=EMA_neg_social_comp-ci, ymax=EMA_neg_social_comp+ci), width=.2,
                 position=position_dodge(0.05))
p + theme_classic() + scale_color_manual(values=c('orange1','olivedrab3'))

# Does it make sense to have an N? Do I use N for Pps or for observations

Modeling approach using a dummy coding - Negative Social Interactions

Visualize cell means and marginal means

library(apaTables)

soc_neg_means <- apa.2way.table(iv1 = Condition, 
               iv2 = onewkpost_14wkpost, 
               dv = EMA_neg_social_comp, 
               data = socialconnect_scale_include2,
               show.marginal.means = TRUE)

soc_neg_means
## 
## 
## Means and standard deviations for EMA_neg_social_comp as a function of a 2(Condition) X 2(onewkpost_14wkpost) design 
## 
##            onewkpost_14wkpost                               
##                         1week      14week      Marginal     
##  Condition                  M   SD      M   SD        M   SD
##          0               1.96 1.22   1.87 1.18     1.91 1.20
##          1               1.93 1.10   1.69 1.13     1.81 1.12
##   Marginal               1.94 1.16   1.78 1.16              
## 
## Note. M and SD represent mean and standard deviation, respectively. 
## Marginal indicates the means and standard deviations pertaining to main effects.

Creating dummy codes to use in the model (See reference note 1.a).

#Dummy code using if/else statement (for Condition variable): 
socialconnect_scale_include2$Dummy_condition <- if_else(socialconnect_scale_include2$Condition == 1, 1, 0) #Dummy code Condition 1 as "1" and Condition 0 as "0"

#Dummy code using if/else statement (for Time variable):
socialconnect_scale_include2$Dummy_week <- if_else(socialconnect_scale_include2$onewkpost_14wkpost == "14week", 1, 0) #Dummy code 1-week post as "0" and 14-weeks post as "1"

#Checking 
#table(socialconnect_scale_include2$Dummy_condition,socialconnect_scale_include2$Dummy_week)

#Running the multi-level model using dummy-coded variables (see Reference note 1.c).

# Modeling for 1-week post 
model_neg_soc_Dummy_1wkPost <- lmer (EMA_neg_social_comp ~  Dummy_condition*Dummy_week +  (1 | PID), #(1 | PID) accounts for dependency
             data = socialconnect_scale_include2)
summary(model_neg_soc_Dummy_1wkPost) #Gives same output as if Condition were not dummy coded.  
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EMA_neg_social_comp ~ Dummy_condition * Dummy_week + (1 | PID)
##    Data: socialconnect_scale_include2
## 
## REML criterion at convergence: 3980.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6894 -0.5969 -0.2590  0.4411  5.0227 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PID      (Intercept) 0.3744   0.6119  
##  Residual             0.9734   0.9866  
## Number of obs: 1345, groups:  PID, 111
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                 1.946e+00  9.941e-02  1.533e+02  19.571   <2e-16
## Dummy_condition             1.502e-03  1.395e-01  1.509e+02   0.011    0.991
## Dummy_week                 -6.995e-02  7.692e-02  1.240e+03  -0.909    0.363
## Dummy_condition:Dummy_week -1.613e-01  1.086e-01  1.244e+03  -1.485    0.138
##                               
## (Intercept)                ***
## Dummy_condition               
## Dummy_week                    
## Dummy_condition:Dummy_week    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dmmy_c Dmmy_w
## Dummy_cndtn -0.713              
## Dummy_week  -0.400  0.285       
## Dmmy_cnd:D_  0.283 -0.390 -0.708
library(sjPlot)
sjPlot::tab_model(model_neg_soc_Dummy_1wkPost)
  EMA neg social comp
Predictors Estimates CI p
(Intercept) 1.95 1.75 – 2.14 <0.001
Dummy_condition 0.00 -0.27 – 0.27 0.991
Dummy_week -0.07 -0.22 – 0.08 0.363
Dummy_condition *
Dummy_week 		-0.16 -0.37 – 0.05 0.137
Random Effects
σ2 0.97
τ00 PID 0.37
ICC 0.28
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.006 / 0.282 
plot_model(model_neg_soc_Dummy_1wkPost, type = "pred", terms = c("Dummy_week[1,0]", "Dummy_condition"))

Using dummy codes to test the null hypotheses that there is no significant difference between C0 and C1 for social connection at 1-week post

fit_neg <- lm(model_neg_soc_Dummy_1wkPost) #Looking at the lm output (instead of anova(model)) for it (because ANOVAs are just fancy linear models) 
summary(fit_neg)
## 
## Call:
## lm(formula = model_neg_soc_Dummy_1wkPost)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1104 -0.8772 -0.6217  0.8915  6.2639 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.7152386  0.1076366  15.935  < 2e-16 ***
## Dummy_condition            -0.0179207  0.0894439  -0.200  0.84123    
## Dummy_week                 -0.0809203  0.0895513  -0.904  0.36636    
## PID                         0.0004849  0.0001738   2.790  0.00535 ** 
## Dummy_condition:Dummy_week -0.1532498  0.1261135  -1.215  0.22451    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.156 on 1340 degrees of freedom
## Multiple R-squared:  0.01329,    Adjusted R-squared:  0.01035 
## F-statistic: 4.514 on 4 and 1340 DF,  p-value: 0.001264
lm(model_neg_soc_Dummy_1wkPost) %>% 
  broom::tidy()
## # A tibble: 5 x 5
##   term                        estimate std.error statistic  p.value
##   <chr>                          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)                 1.72      0.108       15.9   1.74e-52
## 2 Dummy_condition            -0.0179    0.0894      -0.200 8.41e- 1
## 3 Dummy_week                 -0.0809    0.0896      -0.904 3.66e- 1
## 4 PID                         0.000485  0.000174     2.79  5.35e- 3
## 5 Dummy_condition:Dummy_week -0.153     0.126       -1.22  2.25e- 1

#Interpreting the output B0 (intercept) + B1 (Dummy_condition) * onewkpost_14wkpost (Dummy_week) Intercept (B0): The mean for Condition 0 at 1-week post (i.e. when both dummys are zero)
Dummy_week (B1): The difference across time for Condition 0 Dummy_condition (simple effect - the main effect of condition): Codes for the difference between the groups at 1-week post. Here, p = 0.06585

Using dummy codes to test the null hypotheses that there is no significant difference between C0 and C1 for social connection at 14-weeks post

#OPTION

#Dummy code using if/else statement (for Time variable):
socialconnect_scale_include2$Dummy_week <- if_else(socialconnect_scale_include2$onewkpost_14wkpost == "1week", 1, 0) #Dummy code 1-week post as "1" and 14-weeks post as "0"

#Checking 
#table(socialconnect_scale_include2$Dummy_condition,socialconnect_scale_include2$Dummy_week)

# Modeling for 14-week post 
model_neg_soc_Dummy_14wkPost <- lmer (EMA_neg_social_comp ~  Dummy_condition*Dummy_week +  (1 | PID), #(1 | PID) accounts for dependency
             data = socialconnect_scale_include2)
summary(model_neg_soc_Dummy_14wkPost) #Gives same output as if Condition were not dummy coded.  
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EMA_neg_social_comp ~ Dummy_condition * Dummy_week + (1 | PID)
##    Data: socialconnect_scale_include2
## 
## REML criterion at convergence: 3980.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6894 -0.5969 -0.2590  0.4411  5.0227 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  PID      (Intercept) 0.3744   0.6119  
##  Residual             0.9734   0.9866  
## Number of obs: 1345, groups:  PID, 111
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                   1.87557    0.09840  147.43319  19.061   <2e-16
## Dummy_condition              -0.15984    0.13941  150.55508  -1.147    0.253
## Dummy_week                    0.06995    0.07692 1240.12952   0.909    0.363
## Dummy_condition:Dummy_week    0.16134    0.10862 1243.91501   1.485    0.138
##                               
## (Intercept)                ***
## Dummy_condition               
## Dummy_week                    
## Dummy_condition:Dummy_week    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dmmy_c Dmmy_w
## Dummy_cndtn -0.706              
## Dummy_week  -0.378  0.267       
## Dmmy_cnd:D_  0.267 -0.389 -0.708
library(sjPlot)
sjPlot::tab_model(model_neg_soc_Dummy_14wkPost)
  EMA neg social comp
Predictors Estimates CI p
(Intercept) 1.88 1.68 – 2.07 <0.001
Dummy_condition -0.16 -0.43 – 0.11 0.252
Dummy_week 0.07 -0.08 – 0.22 0.363
Dummy_condition *
Dummy_week 		0.16 -0.05 – 0.37 0.137
Random Effects
σ2 0.97
τ00 PID 0.37
ICC 0.28
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.006 / 0.282 
fit_14_neg <- lm(model_neg_soc_Dummy_14wkPost) #Looking at the lm output (instead of anova(model)) for it (because ANOVAs are just fancy linear models) 
summary(fit_14_neg)
## 
## Call:
## lm(formula = model_neg_soc_Dummy_14wkPost)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1104 -0.8772 -0.6217  0.8915  6.2639 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 1.6343183  0.1053328  15.516  < 2e-16 ***
## Dummy_condition            -0.1711705  0.0890097  -1.923  0.05469 .  
## Dummy_week                  0.0809203  0.0895513   0.904  0.36636    
## PID                         0.0004849  0.0001738   2.790  0.00535 ** 
## Dummy_condition:Dummy_week  0.1532498  0.1261135   1.215  0.22451    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.156 on 1340 degrees of freedom
## Multiple R-squared:  0.01329,    Adjusted R-squared:  0.01035 
## F-statistic: 4.514 on 4 and 1340 DF,  p-value: 0.001264
#lm(model_soc_Dummy_14wkPost) %>% 
  #broom::tidy()

MODELS FOR INDIVIDUAL ITEMS

NUMBER OF SOCIAL INTERACTIONS

Describing the data

# Total number of social interactions (soc_numpeople_imp)

# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$soc_numpeople_imp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 323 3.431447 1.140636      3 3.472423 1.4826   1   5     4
## X12    2      1    1 347 3.303624 1.068567      3 3.302357 1.4826   1   5     4
##            skew   kurtosis         se
## X11 -0.21942506 -0.8245707 0.06346670
## X12 -0.05801778 -0.7291598 0.05736367
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$soc_numpeople_imp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 344 3.080047 1.241756      3 3.110639 1.4826   0   5     5
## X12    2      1    1 331 2.903863 1.241911      3 2.902561 1.4826   0   5     5
##            skew   kurtosis         se
## X11 -0.03365631 -0.8264586 0.06695100
## X12 -0.01720361 -0.6383279 0.06826158
# Post
psych::describeBy(socialconnect_scale_include2$soc_numpeople_imp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 667 3.250216 1.205842      3 3.285783 1.4826   0   5     5
## X12    2      1    1 678 3.108461 1.172758      3 3.127824 1.4826   0   5     5
##           skew   kurtosis         se
## X11 -0.1448039 -0.8151606 0.04669038
## X12 -0.1068199 -0.5861024 0.04503951
boxplot(soc_numpeople_imp~Condition,socialconnect_scale_include2, main="Social Interaction by Condition",
   xlab="Condition", ylab="Total number of social interactions")

Mixed model for Number of social interactions (POST)

# NEED TO UPDATE THESE TO REMOVE "onewkpost_14wkpost" FROM (1 + onewkpost_14wkpost | PID) once I confirm with David. 

EQ_MM_Social_model7 <- lmer (soc_numpeople_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) 

sjPlot::tab_model(EQ_MM_Social_model7)
  soc numpeople imp
Predictors Estimates CI p
(Intercept) 3.43 3.21 – 3.66 <0.001
onewkpost_14wkpost
[14week] 		-0.33 -0.52 – -0.14 0.001
Condition -0.09 -0.40 – 0.22 0.570
onewkpost_14wkpost
[14week] * Condition 		-0.05 -0.31 – 0.22 0.722
Random Effects
σ2 0.65
τ00 PID 0.59
τ11 PID.onewkpost_14wkpost14week 0.27
ρ01 PID 0.02
ICC 0.53
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.024 / 0.541 
#This model does not yet account for person-level predictors 

# Plotting model

#library(sjPlot)
#plot_model(model = EQ_MM_Social_model7, type = "int")

# Calculating Cohen's d (effect size) - using method from Jarrett Byrnes described here:     https://stats.stackexchange.com/questions/95054/how-to-get-an-overall-p-value-and-effect-size-for-a-categorical-factor-in-a-mi

# Attempt #1
#r2.corr.mer(EQ_MM_Social_model7)

# Attempt #2
#MuMIn::r.squaredGLMM(EQ_MM_Social_model7)  
# 'R2m' = marginal R squared (proportion of variance explained by the fixed factor(s) alone)
# 'R2c' = conditional R squared (proportion of variance explained by both the fixed and random factors)

# Standardized effect size 
sjstats::std_beta(EQ_MM_Social_model7, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.17573245
## 2           onewkpost_14wkpost14week     -0.29729505
## 3                          Condition     -0.03802438
## 4 onewkpost_14wkpost14week:Condition     -0.02012180

NUMBER OF UNIQUE SOCIAL INTERACTIONS

Describing the data

# Number of unique social interactions (soc_numdifpeople_imp)

# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$soc_numdifpeople_imp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean        sd median  trimmed    mad min max
## X11    1      0    1 323 3.213456 0.9656918      3 3.188981 1.4826   0   5
## X12    2      1    1 347 3.066127 0.9758412      3 3.039234 1.4826   1   5
##     range        skew    kurtosis         se
## X11     5 -0.02211768 -0.01772336 0.05373253
## X12     4  0.09159777 -0.39869228 0.05238590
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$soc_numdifpeople_imp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 344 2.961974 1.091318      3 2.927244 1.4826   0   5     5
## X12    2      1    1 331 2.652487 1.022515      3 2.637634 1.4826   0   5     5
##          skew   kurtosis         se
## X11 0.1288718 -0.2819911 0.05883995
## X12 0.2051521  0.2586399 0.05620252
# Post 
psych::describeBy(socialconnect_scale_include2$soc_numdifpeople_imp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 667 3.083756 1.039263      3 3.053954 1.4826   0   5     5
## X12    2      1    1 678 2.864188 1.019381      3 2.847278 1.4826   0   5     5
##           skew    kurtosis         se
## X11 0.02519652 -0.17335532 0.04024041
## X12 0.11463383 -0.08426436 0.03914910
boxplot(soc_numdifpeople_imp~Condition,socialconnect_scale_include2, main="Brief Social Connection Scores by Condition",
   xlab="Condition", ylab="Total number of unique social interactions")

#THIS LOOKS OFF....

Mixed model for Number of unique social interactions (POST)

# NEED TO UPDATE THESE TO REMOVE "onewkpost_14wkpost" FROM (1 + onewkpost_14wkpost | PID) once I confirm with David. 

EQ_MM_Social_model6 <- lmer (soc_numdifpeople_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) 

sjPlot::tab_model(EQ_MM_Social_model6)
  soc numdifpeople imp
Predictors Estimates CI p
(Intercept) 3.22 3.04 – 3.40 <0.001
onewkpost_14wkpost
[14week] 		-0.25 -0.42 – -0.07 0.005
Condition -0.14 -0.40 – 0.11 0.269
onewkpost_14wkpost
[14week] * Condition 		-0.13 -0.38 – 0.12 0.301
Random Effects
σ2 0.60
τ00 PID 0.37
τ11 PID.onewkpost_14wkpost14week 0.23
ρ01 PID -0.11
ICC 0.43
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.033 / 0.449 
#This model does not yet account for person-level predictors 

# Plotting model

library(sjPlot)
plot_model(model = EQ_MM_Social_model6, type = "int")

# Calculating Cohen's d (effect size) - using method from Jarrett Byrnes described here:     https://stats.stackexchange.com/questions/95054/how-to-get-an-overall-p-value-and-effect-size-for-a-categorical-factor-in-a-mi

# Attempt #1
#r2.corr.mer(EQ_MM_Social_model6)

# Attempt #2
#MuMIn::r.squaredGLMM(EQ_MM_Social_model6)  
# 'R2m' = marginal R squared (proportion of variance explained by the fixed factor(s) alone)
# 'R2c' = conditional R squared (proportion of variance explained by both the fixed and random factors)

# Standardized effect size 
sjstats::std_beta(EQ_MM_Social_model6, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.16728752
## 2           onewkpost_14wkpost14week     -0.30462175
## 3                          Condition     -0.06944586
## 4 onewkpost_14wkpost14week:Condition     -0.06315597
# Between Group Contrasts from the interaction model

#Getting the Simple Means 
est_means_close <- EQ_MM_Social_model6 %>% 
  emmeans::emmeans(specs = c("onewkpost_14wkpost", "Condition"))  

est_means_close #Gives us all four simple means / SEs
##  onewkpost_14wkpost Condition emmean     SE  df lower.CL upper.CL
##  1week                      0   3.22 0.0925 110     3.04     3.40
##  14week                     0   2.97 0.1076 107     2.76     3.18
##  1week                      1   3.07 0.0912 107     2.89     3.26
##  14week                     1   2.69 0.1078 110     2.48     2.91
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
contrast_means_btw_close <- est_means_close %>% 
  emmeans::contrast(by = "onewkpost_14wkpost", simple = "Condition")

contrast_means_btw_close #Gives contrast stats for simple main effect at pre and then at post, with significance values
## onewkpost_14wkpost = 1week:
##  contrast estimate     SE  df t.ratio p.value
##  0 effect   0.0718 0.0650 109  1.106  0.2713 
##  1 effect  -0.0718 0.0650 109 -1.106  0.2713 
## 
## onewkpost_14wkpost = 14week:
##  contrast estimate     SE  df t.ratio p.value
##  0 effect   0.1372 0.0761 108  1.802  0.0744 
##  1 effect  -0.1372 0.0761 108 -1.802  0.0744 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: fdr method for 2 tests

Running models with outliers removed

Detecting Outliers

#Obtaining values of outliers 

library(dplyr)
library(ggplot2)

is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

socialconnect_scale_include2 %>%
  group_by(Condition) %>%
  mutate(outlier = ifelse(is_outlier(soc_numdifpeople_imp), soc_numdifpeople_imp, as.numeric(NA))) %>%
  ggplot(., aes(x = factor(Condition), y = soc_numdifpeople_imp)) +
    geom_boxplot() +
    geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)

# Removing Outliers

socialconnect_out1 <- socialconnect_scale_include2 %>% 
  filter(Condition == 0, soc_numdifpeople_imp > 0) 

#View(socialconnect_scale_include2)

socialconnect_out2 <- socialconnect_scale_include2 %>% 
  filter(Condition == 1, 5 > soc_numdifpeople_imp, soc_numdifpeople_imp > 0) 

socialconnect_out <- rbind(socialconnect_out1, socialconnect_out2)

Re-running Mixed model for Number of unique social interactions (POST, outliers removed)

# NEED TO UPDATE THESE TO REMOVE "onewkpost_14wkpost" FROM (1 + onewkpost_14wkpost | PID) once I confirm with David. 

EQ_MM_Social_model_out <- lmer (soc_numdifpeople_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_out) 

sjPlot::tab_model(EQ_MM_Social_model_out)
  soc numdifpeople imp
Predictors Estimates CI p
(Intercept) 3.23 3.06 – 3.39 <0.001
onewkpost_14wkpost
[14week] 		-0.23 -0.38 – -0.09 0.002
Condition -0.31 -0.54 – -0.08 0.009
onewkpost_14wkpost
[14week] * Condition 		-0.11 -0.31 – 0.10 0.320
Random Effects
σ2 0.51
τ00 PID 0.30
τ11 PID.onewkpost_14wkpost14week 0.13
ρ01 PID -0.07
ICC 0.41
N PID 110
Observations 1290
Marginal R2 / Conditional R2 0.058 / 0.441 
#This model does not yet account for person-level predictors 

# Plotting model

library(sjPlot)
plot_model(model = EQ_MM_Social_model_out, type = "int")

/// SOCIAL CLOSENESS / CONNECTION ///

# Descriptives - Social closeness/connection (soc_closeness_imp)

# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$soc_closeness_imp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 323 2.029290 1.781433      3 2.264327 1.4826  -3   4     7
## X12    2      1    1 347 1.625134 2.027648      3 1.874270 1.4826  -3   4     7
##          skew   kurtosis         se
## X11 -1.075461  0.5779050 0.09912159
## X12 -0.906075 -0.1162136 0.10884983
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$soc_closeness_imp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 344 1.823230 1.970237      3 2.076779 1.4826  -3   4     7
## X12    2      1    1 331 1.452358 2.210033      1 1.682001 2.9652  -3   4     7
##           skew   kurtosis        se
## X11 -0.9601319  0.1325495 0.1062281
## X12 -0.6562069 -0.7061560 0.1214744
# Post
psych::describeBy(socialconnect_scale_include2$soc_closeness_imp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd   median  trimmed      mad min max
## X11    1      0    1 667 1.923016 1.882594 3.000000 2.167573 1.482600  -3   4
## X12    2      1    1 678 1.540784 2.118843 1.730337 1.780610 1.882402  -3   4
##     range       skew   kurtosis         se
## X11     7 -1.0252306  0.3657528 0.07289433
## X12     7 -0.7804077 -0.4308346 0.08137366
boxplot(soc_closeness_imp~Condition,socialconnect_scale_include2, main="Brief Social Connection Scores by Condition",
   xlab="Condition", ylab="Perceived social closeness / connection")

Mixed model for social closeness (POST)

#Model 4 social closeness
EQ_MM_Social_model_close <- lmer (soc_closeness_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) 
 
#this is a group by time interaction
#3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)
# "1 + onewkpost_14wkpost | PID" estimate a random intercept and a random slope for each individual 
summary(EQ_MM_Social_model_close)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## soc_closeness_imp ~ onewkpost_14wkpost * Condition + (1 + onewkpost_14wkpost |  
##     PID)
##    Data: socialconnect_scale_include2
## 
## REML criterion at convergence: 5452.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1960 -0.4915  0.1954  0.6618  3.2262 
## 
## Random effects:
##  Groups   Name                     Variance Std.Dev. Corr
##  PID      (Intercept)              0.8587   0.9267       
##           onewkpost_14wkpost14week 0.5958   0.7719   0.07
##  Residual                          2.8257   1.6810       
## Number of obs: 1345, groups:  PID, 111
## 
## Fixed effects:
##                                     Estimate Std. Error        df t value
## (Intercept)                          2.05936    0.15678 113.45533  13.136
## onewkpost_14wkpost14week            -0.23918    0.16792 104.66484  -1.424
## Condition                           -0.40727    0.21975 111.25272  -1.853
## onewkpost_14wkpost14week:Condition   0.03363    0.23742 104.35512   0.142
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## onewkpost_14wkpost14week             0.1573    
## Condition                            0.0665 .  
## onewkpost_14wkpost14week:Condition   0.8876    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) on_1414 Condtn
## onwkps_1414 -0.305               
## Condition   -0.713  0.217        
## onwk_1414:C  0.216 -0.707  -0.296
sjPlot::tab_model(EQ_MM_Social_model_close)
  soc closeness imp
Predictors Estimates CI p
(Intercept) 2.06 1.75 – 2.37 <0.001
onewkpost_14wkpost
[14week] 		-0.24 -0.57 – 0.09 0.154
Condition -0.41 -0.84 – 0.02 0.064
onewkpost_14wkpost
[14week] * Condition 		0.03 -0.43 – 0.50 0.887
Random Effects
σ2 2.83
τ00 PID 0.86
τ11 PID.onewkpost_14wkpost14week 0.60
ρ01 PID 0.07
ICC 0.30
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.012 / 0.308 
#This model does not yet account for person-level predictors

#Model 4 
library(sjPlot)
plot_model(model = EQ_MM_Social_model_close, type = "int")

#This model does not yet account for person-level predictors 

# Calculating Cohen's d (effect size) - using method from Jarrett Byrnes described here:     https://stats.stackexchange.com/questions/95054/how-to-get-an-overall-p-value-and-effect-size-for-a-categorical-factor-in-a-mi

# Attempt #1
#r2.corr.mer(EQ_MM_Social_model_close)

# Attempt #2
#MuMIn::r.squaredGLMM(EQ_MM_Social_model_close)  
# 'R2m' = marginal R squared (proportion of variance explained by the fixed factor(s) alone)
# 'R2c' = conditional R squared (proportion of variance explained by both the fixed and random factors)

# Standardized effect size 
sjstats::std_beta(EQ_MM_Social_model_close, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)     0.061444628
## 2           onewkpost_14wkpost14week    -0.110367333
## 3                          Condition    -0.101168294
## 4 onewkpost_14wkpost14week:Condition     0.008354016

Removing outliers

#Obtaining values of outliers 

library(dplyr)
library(ggplot2)

is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

socialconnect_scale_include2 %>%
  group_by(Condition) %>%
  mutate(outlier = ifelse(is_outlier(soc_closeness_imp), soc_closeness_imp, as.numeric(NA))) %>%
  ggplot(., aes(x = factor(Condition), y = soc_closeness_imp)) +
    geom_boxplot() +
    geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)

# NEED TO FIGURE OUT HOW TO Remove outliers at both levels

Removing Outliers

socialconnect_scale_include2_out4 <- socialconnect_scale_include2 %>% 
  filter(Condition == 0, soc_closeness_imp > -3) 

socialconnect_scale_include2_1_4 <- socialconnect_scale_include2 %>% 
  filter(Condition == 1, soc_closeness_imp > -3) 

socialconnect_scale_out4 <- rbind(socialconnect_scale_include2_out4, socialconnect_scale_include2_1_4)
#run the models with and without these outliers

Mixed model for EMA Social Connection - Closeness (outliers removed)

library(lme4)
library(lmerTest)

socialconnect_scale_out4 <- socialconnect_scale_out4 %>% 
  mutate(onewkpost_14wkpost = as.factor(onewkpost_14wkpost), 
         timepoint = as.factor(timepoint))

EQ_MM_SocialConnect_model_out4 <- lmer (soc_closeness_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_out4) 
 
#this is a group by time interaction
#3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)
# "1 + onewkpost_14wkpost | PID" estimate a random intercept and a random slope for each individual 
#summary(EQ_MM_SocialConnect_model_out4)
sjPlot::tab_model(EQ_MM_SocialConnect_model_out4)
  soc closeness imp
Predictors Estimates CI p
(Intercept) 2.22 1.99 – 2.45 <0.001
onewkpost_14wkpost
[14week] 		-0.12 -0.41 – 0.18 0.440
Condition -0.22 -0.55 – 0.11 0.182
onewkpost_14wkpost
[14week] * Condition 		-0.05 -0.46 – 0.37 0.833
Random Effects
σ2 1.90
τ00 PID 0.44
τ11 PID.onewkpost_14wkpost14week 0.55
ρ01 PID 0.04
ICC 0.28
N PID 111
Observations 1253
Marginal R2 / Conditional R2 0.007 / 0.283 
#This model does not yet account for person-level predictors 

# Standardized effect size 
sjstats::std_beta(EQ_MM_SocialConnect_model_out4, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.01932192
## 2           onewkpost_14wkpost14week     -0.08585512
## 3                          Condition     -0.06974864
## 4 onewkpost_14wkpost14week:Condition     -0.01402535
library(sjPlot)
plot_model(model = EQ_MM_SocialConnect_model_out4, type = "int")

Creating separate dataframes for 1-week and 14-week follow-up (outliers removed)

socialconnect_scale_out_1wkpost <- socialconnect_scale_out4 %>% 
  filter(onewkpost_14wkpost == '1week')

socialconnect_scale_out_14wkpost <- socialconnect_scale_out4 %>% 
  filter(onewkpost_14wkpost == '14week')

###Check to see if outliers were adequately removed.

/// POSITIVE SOCIAL INTERACTIONS ///

# Number of positive social interactions (soc_positive_imp) 

# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$soc_positive_imp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 323 2.697095 1.116944      3 2.718771 1.4826   0   5     5
## X12    2      1    1 347 2.666172 1.183559      3 2.681583 1.4826   0   5     5
##            skew   kurtosis         se
## X11 -0.06471052 -0.5336181 0.06214843
## X12 -0.06403106 -0.6010315 0.06353680
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$soc_positive_imp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 344 2.553266 1.277810      2 2.475809 1.4826   0   5     5
## X12    2      1    1 331 2.406589 1.338948      2 2.398418 1.4826   0   5     5
##           skew   kurtosis         se
## X11 0.34129538 -0.5951527 0.06889492
## X12 0.07087303 -0.7941901 0.07359526
# Post
psych::describeBy(socialconnect_scale_include2$soc_positive_imp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n     mean       sd median  trimmed    mad min max range
## X11    1      0    1 667 2.622916 1.203855      3 2.591561 1.4826   0   5     5
## X12    2      1    1 678 2.539443 1.267544      3 2.543644 1.4826   0   5     5
##            skew   kurtosis         se
## X11  0.16173932 -0.5744157 0.04661346
## X12 -0.02326539 -0.6910354 0.04867975
boxplot(soc_positive_imp~Condition,socialconnect_scale_include2, main="Brief Social Connection Scores by Condition",
   xlab="Condition", ylab="Positive Interactions")

Mixed model for positive social interactions (POST)

#Model positive social interactions
EQ_MM_Social_model_pos <- lmer (soc_positive_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) 
 
#this is a group by time interaction
#3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)
# "1 + onewkpost_14wkpost | PID" estimate a random intercept and a random slope for each individual 
#summary(EQ_MM_Social_model_close)
sjPlot::tab_model(EQ_MM_Social_model_pos)
  soc positive imp
Predictors Estimates CI p
(Intercept) 2.70 2.46 – 2.94 <0.001
onewkpost_14wkpost
[14week] 		-0.14 -0.31 – 0.03 0.112
Condition -0.02 -0.36 – 0.32 0.901
onewkpost_14wkpost
[14week] * Condition 		-0.06 -0.31 – 0.18 0.600
Random Effects
σ2 0.64
τ00 PID 0.71
τ11 PID.onewkpost_14wkpost14week 0.20
ρ01 PID 0.18
ICC 0.58
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.005 / 0.581 
#This model does not yet account for person-level predictors

#Model 4 
library(sjPlot)
plot_model(model = EQ_MM_Social_model_pos, type = "int")

# Standardized effect size 
sjstats::std_beta(EQ_MM_Social_model_pos, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)     0.087737856
## 2           onewkpost_14wkpost14week    -0.138550060
## 3                          Condition    -0.008669619
## 4 onewkpost_14wkpost14week:Condition    -0.026209135

Removing Outliers

/// DEGREE OF POSITIVE SOCIAL INTERACTIONS ///

# Degree of positive social interactions (soc_howpos_imp)

# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$soc_howpos_imp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean        sd  median  trimmed       mad min max
## X11    1      0    1 323 2.200946 0.4747391 2.21063 2.200407 0.3122799   1   3
## X12    2      1    1 347 2.241808 0.5133123 2.21063 2.264901 0.3122799   1   3
##     range        skew  kurtosis         se
## X11     2 -0.06938184 0.9631913 0.02641519
## X12     2 -0.25338194 0.5105925 0.02755605
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$soc_howpos_imp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean        sd median  trimmed       mad min max
## X11    1      0    1 344 2.231654 0.5278829      2 2.241627 0.0000000   1   3
## X12    2      1    1 331 2.165545 0.5784341      2 2.206775 0.3122799   1   3
##     range        skew    kurtosis         se
## X11     2 -0.00536284  0.01646169 0.02846154
## X12     2 -0.20002055 -0.10722114 0.03179361
# Post
psych::describeBy(socialconnect_scale_include2$soc_howpos_imp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n     mean        sd  median  trimmed       mad min max
## X11    1      0    1 667 2.216783 0.5027090 2.00000 2.221672 0.3122799   1   3
## X12    2      1    1 678 2.204576 0.5470021 2.21063 2.243939 0.3122799   1   3
##     range        skew  kurtosis         se
## X11     2 -0.02187152 0.4266876 0.01946497
## X12     2 -0.24739523 0.1921725 0.02100749
boxplot(soc_howpos_imp~Condition,socialconnect_scale_include2, main="Brief Social Connection Scores by Condition",
   xlab="Condition", ylab="Degree of Positive Interactions")

Mixed model for Degree of positive social interactions (POST)

#Model degree of positive social interactions
EQ_MM_Social_model9 <- lmer (soc_howpos_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) 
 
sjPlot::tab_model(EQ_MM_Social_model9)
  soc howpos imp
Predictors Estimates CI p
(Intercept) 2.21 2.14 – 2.28 <0.001
onewkpost_14wkpost
[14week] 		0.03 -0.06 – 0.11 0.535
Condition 0.04 -0.06 – 0.14 0.438
onewkpost_14wkpost
[14week] * Condition 		-0.10 -0.22 – 0.02 0.097
Random Effects
σ2 0.22
τ00 PID 0.03
τ11 PID.onewkpost_14wkpost14week 0.03
ρ01 PID 0.33
ICC 0.21
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.003 / 0.213 
#This model does not yet account for person-level predictors

#Model 4 
library(sjPlot)
plot_model(model = EQ_MM_Social_model9, type = "int")

# Standardized effect size 
sjstats::std_beta(EQ_MM_Social_model9, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.03043853
## 2           onewkpost_14wkpost14week     -0.04586647
## 3                          Condition      0.03740360
## 4 onewkpost_14wkpost14week:Condition     -0.09560889

Removing outliers

#Obtaining values of outliers 

library(dplyr)
library(ggplot2)

is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

socialconnect_scale_include2 %>%
  group_by(Condition) %>%
  mutate(outlier = ifelse(is_outlier(soc_howpos_imp), soc_howpos_imp, as.numeric(NA))) %>%
  ggplot(., aes(x = factor(Condition), y = soc_howpos_imp)) +
    geom_boxplot() +
    geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)

# NEED TO FIGURE OUT HOW TO Remove outliers at both levels

Removing Outliers

socialconnect_scale_include2_out5 <- socialconnect_scale_include2 %>% 
  filter(Condition == 0, soc_howpos_imp > 1) %>% 
  filter(Condition == 0, soc_howpos_imp < 3) 

socialconnect_scale_include2_1_5 <- socialconnect_scale_include2 %>% 
  filter(Condition == 1, soc_howpos_imp > 1) %>% 
  filter(Condition == 1, soc_howpos_imp < 3)

socialconnect_scale_out5 <- rbind(socialconnect_scale_include2_out5, socialconnect_scale_include2_1_5)
#run the models with and without these outliers

Mixed model for EMA Social Connection - Closeness (outliers removed)

library(lme4)
library(lmerTest)

socialconnect_scale_out5 <- socialconnect_scale_out4 %>% 
  mutate(onewkpost_14wkpost = as.factor(onewkpost_14wkpost), 
         timepoint = as.factor(timepoint))

EQ_MM_SocialConnect_model_out5 <- lmer (soc_howpos_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_out5) 
 
#this is a group by time interaction
#3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)
# "1 + onewkpost_14wkpost | PID" estimate a random intercept and a random slope for each individual 
#summary(EQ_MM_SocialConnect_model_out4)
sjPlot::tab_model(EQ_MM_SocialConnect_model_out5)
  soc howpos imp
Predictors Estimates CI p
(Intercept) 2.23 2.18 – 2.28 <0.001
onewkpost_14wkpost
[14week] 		0.03 -0.07 – 0.13 0.601
Condition 0.05 -0.02 – 0.12 0.188
onewkpost_14wkpost
[14week] * Condition 		-0.09 -0.23 – 0.05 0.222
Random Effects
σ2 0.22
τ00 PID 0.00
τ11 PID.onewkpost_14wkpost14week 0.07
ρ01 PID  
N PID 111
Observations 1253
Marginal R2 / Conditional R2 0.003 / NA 
#This model does not yet account for person-level predictors 

# Standardized effect size 
sjstats::std_beta(EQ_MM_SocialConnect_model_out5, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.01374715
## 2           onewkpost_14wkpost14week     -0.03462212
## 3                          Condition      0.04880545
## 4 onewkpost_14wkpost14week:Condition     -0.08948873
library(sjPlot)
plot_model(model = EQ_MM_SocialConnect_model_out5, type = "int")

/// NEGATIVE SOCIAL INTERACTIONS ///

# Number of negative social interactions (soc_neg_imp)

# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$soc_neg_imp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n      mean        sd median   trimmed mad min max range
## X11    1      0    1 323 0.7065164 0.9232930      0 0.5374703   0   0   5     5
## X12    2      1    1 347 0.6219228 0.7769684      0 0.4939326   0   0   4     4
##         skew kurtosis         se
## X11 1.556559 2.616528 0.05137340
## X12 1.213313 1.258263 0.04170985
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$soc_neg_imp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n      mean        sd median   trimmed mad min max range
## X11    1      0    1 344 0.6174698 0.8644165      0 0.4616290   0   0   5     5
## X12    2      1    1 331 0.4646743 0.7969989      0 0.3011593   0   0   5     5
##         skew kurtosis         se
## X11 1.580285 2.805182 0.04660621
## X12 2.477937 8.423070 0.04380702
# Post
psych::describeBy(socialconnect_scale_include2$soc_neg_imp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n      mean        sd median   trimmed mad min max range
## X11    1      0    1 667 0.6605913 0.8938489      0 0.4983447   0   0   5     5
## X12    2      1    1 678 0.5451540 0.7901541      0 0.3889971   0   0   5     5
##         skew kurtosis         se
## X11 1.577675 2.758382 0.03460997
## X12 1.824306 4.592551 0.03034569
boxplot(soc_neg_imp~Condition,socialconnect_scale_include2, main="Brief Social Connection Scores by Condition",
   xlab="Condition", ylab="Negative Interactions")

#WHY ARE THESE SO SIMILAR?

Mixed model for negative social interactions (POST)

#Model negative social interactions
EQ_MM_Social_model_neg <- lmer (soc_neg_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) 
 
#this is a group by time interaction
#3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)
# "1 + onewkpost_14wkpost | PID" estimate a random intercept and a random slope for each individual 
#summary(EQ_MM_Social_model_close)
sjPlot::tab_model(EQ_MM_Social_model_neg)
  soc neg imp
Predictors Estimates CI p
(Intercept) 0.70 0.55 – 0.85 <0.001
onewkpost_14wkpost
[14week] 		-0.08 -0.24 – 0.08 0.320
Condition -0.06 -0.27 – 0.15 0.585
onewkpost_14wkpost
[14week] * Condition 		-0.07 -0.29 – 0.15 0.556
Random Effects
σ2 0.47
τ00 PID 0.25
τ11 PID.onewkpost_14wkpost14week 0.19
ρ01 PID -0.43
ICC 0.35
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.008 / 0.350 
#This model does not yet account for person-level predictors

#Model 4 
library(sjPlot)
plot_model(model = EQ_MM_Social_model_neg, type = "int")

# Standardized effect size 
sjstats::std_beta(EQ_MM_Social_model_neg, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.07884265
## 2           onewkpost_14wkpost14week     -0.13365482
## 3                          Condition     -0.03517613
## 4 onewkpost_14wkpost14week:Condition     -0.03933635

Regression for negative social interactions (1-WEEK POST)

EQ_MM_Social_model5_1wk <- lm (soc_neg_imp ~ Condition, data = socialconnect_scale_1wkpost) 

library(sjPlot)
sjPlot::tab_model(EQ_MM_Social_model5_1wk)
  soc neg imp
Predictors Estimates CI p
(Intercept) 0.71 0.61 – 0.80 <0.001
Condition -0.08 -0.21 – 0.04 0.199
Observations 670
R2 / R2 adjusted 0.002 / 0.001

Regression for negative social interactions (14-WEEK POST)

EQ_MM_Social_model5_14wk <- lm (soc_neg_imp ~ Condition, data = socialconnect_scale_14wkpost) 

library(sjPlot)
sjPlot::tab_model(EQ_MM_Social_model5_14wk)
  soc neg imp
Predictors Estimates CI p
(Intercept) 0.62 0.53 – 0.71 <0.001
Condition -0.15 -0.28 – -0.03 0.017
Observations 675
R2 / R2 adjusted 0.008 / 0.007

Removing outliers

/// DEGREE OF NEGATIVE SOCIAL INTERACTIONS ///

# Degree of negative social interactions (soc_howneg_imp)

# 1-week post
psych::describeBy(socialconnect_scale_1wkpost$soc_howneg_imp, socialconnect_scale_1wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean        sd median  trimmed mad min max range
## X11    1      0    1 323 1.249532 0.5063675      1 1.141308   0   1   3     2
## X12    2      1    1 347 1.304954 0.5327504      1 1.210821   0   1   3     2
##         skew kurtosis         se
## X11 1.987501 3.157264 0.02817504
## X12 1.571219 1.567486 0.02859954
# 14-weeks post
psych::describeBy(socialconnect_scale_14wkpost$soc_howneg_imp, socialconnect_scale_14wkpost$Condition, mat = TRUE)
##     item group1 vars   n     mean        sd median  trimmed mad min max range
## X11    1      0    1 344 1.253777 0.5127616      1 1.146012   0   1   3     2
## X12    2      1    1 331 1.229061 0.4957564      1 1.116299   0   1   3     2
##         skew kurtosis         se
## X11 1.936428 2.911768 0.02764625
## X12 2.142560 3.822752 0.02724924
# Post
psych::describeBy(socialconnect_scale_include2$soc_howneg_imp, socialconnect_scale_include2$Condition, mat = TRUE)
##     item group1 vars   n     mean        sd median  trimmed mad min max range
## X11    1      0    1 667 1.251721 0.5092972      1 1.143735   0   1   3     2
## X12    2      1    1 678 1.267903 0.5160411      1 1.164776   0   1   3     2
##         skew kurtosis         se
## X11 1.965323 3.046739 0.01972007
## X12 1.827956 2.499665 0.01981844
boxplot(soc_howneg_imp~Condition,socialconnect_scale_include2, main="Brief Social Connection Scores by Condition",
   xlab="Condition", ylab="Degree of Negative Interactions")

#WHY ARE THESE LIKE THIS?

Mixed model for degree of negative social interactions (POST)

EQ_MM_Social_model_negd <- lmer (soc_howneg_imp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_include2) 

sjPlot::tab_model(EQ_MM_Social_model_negd)
  soc howneg imp
Predictors Estimates CI p
(Intercept) 1.25 1.17 – 1.32 <0.001
onewkpost_14wkpost
[14week] 		0.01 -0.07 – 0.08 0.865
Condition 0.06 -0.04 – 0.16 0.253
onewkpost_14wkpost
[14week] * Condition 		-0.08 -0.19 – 0.02 0.126
Random Effects
σ2 0.22
τ00 PID 0.04
τ11 PID.onewkpost_14wkpost14week 0.01
ρ01 PID -0.17
ICC 0.16
N PID 111
Observations 1345
Marginal R2 / Conditional R2 0.003 / 0.160 
#This model does not yet account for person-level predictors

#Model 4 
library(sjPlot)
plot_model(model = EQ_MM_Social_model_negd, type = "int")

# Standardized effect size 
sjstats::std_beta(EQ_MM_Social_model_negd, ci.lvl = 0.95,) #"std.estimate" represents the standardized beta coefficients.
##                            Parameter Std_Coefficient
## 1                        (Intercept)      0.03684453
## 2           onewkpost_14wkpost14week     -0.06951643
## 3                          Condition      0.05901677
## 4 onewkpost_14wkpost14week:Condition     -0.08175094

Removing Outliers

******BELOW NEEDS TO BE UPDATED AND ADDED TO NEG COMP SECTION - CODE IS FROM COMP 1****** # Removing Outliers

#Obtaining values of outliers 

library(dplyr)
library(ggplot2)

is_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

socialconnect_scale_include2 %>%
  group_by(Condition) %>%
  mutate(outlier = ifelse(is_outlier(EMA_social_comp), EMA_social_comp, as.numeric(NA))) %>%
  ggplot(., aes(x = factor(Condition), y = EMA_social_comp)) +
    geom_boxplot() +
    geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)

# NEED TO FIGURE OUT HOW TO Remove outliers at both levels

Removing Outliers

socialconnect_scale_include2_out <- socialconnect_scale_include2 %>% 
  filter(Condition == 0, EMA_social_comp > 3) 

socialconnect_scale_include2_1 <- socialconnect_scale_include2 %>% 
  filter(Condition == 1, EMA_social_comp > -2) 

socialconnect_scale_out <- rbind(socialconnect_scale_include2_out, socialconnect_scale_include2_1)
#run the models with and without these outliers

Mixed model for EMA Social Connection (outliers removed)

library(lme4)
library(lmerTest)

socialconnect_scale_out <- socialconnect_scale_out %>% 
  mutate(onewkpost_14wkpost = as.factor(onewkpost_14wkpost), 
         timepoint = as.factor(timepoint))

EQ_MM_SocialConnect_model_out <- lmer (EMA_social_comp ~ onewkpost_14wkpost * Condition + #fixed effects
      (1 + onewkpost_14wkpost | PID), data = socialconnect_scale_out) 
 
#this is a group by time interaction
#3 predictors: 1) Condition(Group), 2) Time, 3) time by condition (group x time)
# "1 + onewkpost_14wkpost | PID" estimate a random intercept and a random slope for each individual 
#summary(EQ_MM_SocialConnect_model_out)
sjPlot::tab_model(EQ_MM_SocialConnect_model_out)
  EMA social comp
Predictors Estimates CI p
(Intercept) 13.68 12.90 – 14.47 <0.001
onewkpost_14wkpost
[14week] 		-0.73 -1.33 – -0.12 0.019
Condition -0.67 -1.77 – 0.43 0.230
onewkpost_14wkpost
[14week] * Condition 		-0.45 -1.31 – 0.41 0.305
Random Effects
σ2 9.74
τ00 PID 7.06
τ11 PID.onewkpost_14wkpost14week 1.92
ρ01 PID 0.29
ICC 0.48
N PID 111
Observations 1328
Marginal R2 / Conditional R2 0.022 / 0.494 
#This model does not yet account for person-level predictors 
library(sjPlot)
plot_model(model = EQ_MM_SocialConnect_model_out, type = "int")

Creating separate dataframes for 1-week and 14-week follow-up (outliers removed)

socialconnect_scale_out_1wkpost <- socialconnect_scale_out %>% 
  filter(onewkpost_14wkpost == '1week')

socialconnect_scale_out_14wkpost <- socialconnect_scale_out %>% 
  filter(onewkpost_14wkpost == '14week')

Regression for Social Connection Composite (1-WEEK POST / outliers removed)

#w/ outliers includeed 
EQ_MM_Social_model_1wk <- lm (EMA_social_comp ~ Condition, data = socialconnect_scale_1wkpost) 
summary(EQ_MM_Social_model_1wk)
## 
## Call:
## lm(formula = EMA_social_comp ~ Condition, data = socialconnect_scale_1wkpost)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.5722  -2.5722   0.3678   3.3078   9.0971 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.5722     0.2308  58.793   <2e-16 ***
## Condition    -0.6694     0.3208  -2.087   0.0373 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.149 on 668 degrees of freedom
## Multiple R-squared:  0.006477,   Adjusted R-squared:  0.004989 
## F-statistic: 4.355 on 1 and 668 DF,  p-value: 0.03729
library(sjPlot)
sjPlot::tab_model(EQ_MM_Social_model_1wk)
  EMA social comp
Predictors Estimates CI p
(Intercept) 13.57 13.12 – 14.03 <0.001
Condition -0.67 -1.30 – -0.04 0.037
Observations 670
R2 / R2 adjusted 0.006 / 0.005 
#w/ outliers removed  
EQ_MM_Social_model_out_1wk <- lm (EMA_social_comp ~ Condition, data = socialconnect_scale_out_1wkpost) 

library(sjPlot)
sjPlot::tab_model(EQ_MM_Social_model_out_1wk)
  EMA social comp
Predictors Estimates CI p
(Intercept) 13.65 13.21 – 14.10 <0.001
Condition -0.75 -1.37 – -0.13 0.018
Observations 668
R2 / R2 adjusted 0.008 / 0.007

Regression for Social Connection Composite (14-WEEK POST / outliers removed) (Note: unclear whether this is a meaningful composite score)

EQ_MM_Social_model_14wk <- lm (EMA_social_comp ~ Condition, data = socialconnect_scale_out_14wkpost) 

library(sjPlot)
sjPlot::tab_model(EQ_MM_Social_model_14wk)
  EMA social comp
Predictors Estimates CI p
(Intercept) 13.06 12.57 – 13.54 <0.001
Condition -1.35 -2.04 – -0.67 <0.001
Observations 660
R2 / R2 adjusted 0.022 / 0.021