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)

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

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

Negative Interactions

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

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

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 == 'Base')

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

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))

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)

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_14wkpostEnd 1.511    0.37
##  Residual                       3.229        
## Number of obs: 1345, groups:  PID, 111
## Fixed Effects:
##                     (Intercept)            onewkpost_14wkpostEnd  
##                         13.6231                          -0.9319  
##                       Condition  onewkpost_14wkpostEnd: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_14wkpostEnd  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 Pr(>|t|)
## (Intercept)                      13.6231     0.4060 110.2270  33.551  < 2e-16
## onewkpost_14wkpostEnd            -0.9319     0.3251 103.6763  -2.867  0.00502
## Condition                        -0.6194     0.5704 109.1143  -1.086  0.27990
## onewkpost_14wkpostEnd:Condition  -0.3074     0.4598 103.2911  -0.669  0.50522
##                                    
## (Intercept)                     ***
## onewkpost_14wkpostEnd           ** 
## Condition                          
## onewkpost_14wkpostEnd:Condition    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) on_14E Condtn
## onwkpst_14E -0.045              
## Condition   -0.712  0.032       
## onwkp_14E: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 [End]	-0.93	-1.57 – -0.29	0.004
Condition	-0.62	-1.74 – 0.50	0.278
onewkpost_14wkpost [End] * Condition	-0.31	-1.21 – 0.59	0.504
Random Effects
σ²	10.43
τ₀₀ _PID	7.24
τ₁₁ _{PID.onewkpost_14wkpostEnd}	2.28
ρ₀₁ _PID	0.37
ICC	0.49
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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).

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

# STATA MLM model EXAMPLE: "mixed outcome_name i.timepoint##i.condition || subjectid:" 
# mixed is the MLM command
# outcome; ## specifies main effects for time and condition as well as interaction effects; 
# || subjectid: is the grouping variable to show that multiple observations are nested within subject
# STATA stats for main effects and interactions: "contrast i.condition##i.timepoint"

# NEED TO ADD code to get contrast means/SEs/stats at each time point (e.g., does condition 1 differ from condition 2 at time 2, accounting for time 1 levels?). @bren: breaking up main effect into simple effects 
## STATA: "margins r.condition#r.timepoint"
# NEED TO ADD code to get within-group means/SEs for change over time
## STATA: "margins R.timepoint, over(condition)"
# NEED TO ADD code to lists means/SEs at each time point
## STATA: "margins timepoint, by(condition)"

#### 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.

# Between Group Contrasts from the interaction model

#install.packages("emmeans")

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

est_means #Gives us all four simple means / SEs

##  onewkpost_14wkpost Condition emmean    SE  df lower.CL upper.CL
##  Base                       0   13.6 0.406 110     12.8     14.4
##  End                        0   12.7 0.509 107     11.7     13.7
##  Base                       1   13.0 0.401 108     12.2     13.8
##  End                        1   11.8 0.508 110     10.8     12.8
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

contrast_means_btw <- est_means %>% 
  emmeans::contrast(by = "onewkpost_14wkpost", simple = "Condition") # Calculate simple effect of condition (aka the difference between Condition 1 and Condition 2 at 1-week post AND the difference between Condition 1 and Condition 2 at 14-weeks post)

contrast_means_btw #Gives contrast stats for simple main effect at pre and then at post, with significance values

## onewkpost_14wkpost = Base:
##  contrast estimate    SE  df t.ratio p.value
##  0 effect    0.310 0.285 109  1.086  0.2800 
##  1 effect   -0.310 0.285 109 -1.086  0.2800 
## 
## onewkpost_14wkpost = End:
##  contrast estimate    SE  df t.ratio p.value
##  0 effect    0.463 0.359 109  1.289  0.2001 
##  1 effect   -0.463 0.359 109 -1.289  0.2001 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: fdr method for 2 tests

Interpreting the simple means (outputs from est_means)

Group 1 1-week post (M = 13.6, SE = 0.41) & end of semester (M = 12.7, SE = 0.51) Group 2 1-week post (M = 13.0, SE = 0.40) & end of semester (M = 11.8, SE = 0.51)

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

The contrast score (difference in means) between group 1 and group 2 was not significantly different at either 1-week post (M = -0.310, SE = 0.28, p = 0.28) or 14-weeks post (M = -0.46, SE = 0.36, p = 0.20).

#QUESTION: Why doesn’t this match up with the reported means? Cameron M thinks there is a good reason for this (bc it is adjusting for the covariate).

Remaining questions - confirm how is this method dealing with the nested structure.

# Within Group Contrasts from the interaction model 
contrast_means_wthn <- est_means %>% 
  emmeans::contrast(by = "Condition", simple = "onewkpost_14wkpost") # Calculate simple effect of time (aka the difference between 1-week post and 14-weeks post for Condition 1 AND the difference between 1-week post and 14-weeks post for Condition 2)

contrast_means_wthn

## Condition = 0:
##  contrast    estimate    SE  df t.ratio p.value
##  Base effect    0.466 0.163 107  2.865  0.0050 
##  End effect    -0.466 0.163 107 -2.865  0.0050 
## 
## Condition = 1:
##  contrast    estimate    SE  df t.ratio p.value
##  Base effect    0.620 0.163 106  3.809  0.0002 
##  End effect    -0.620 0.163 106 -3.809  0.0002 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: fdr method for 2 tests

#For each group separately, there is a significant difference from pre to post.

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

For condition 1, the constrast score (difference in means) between 1-week post and 14-week post was significantly different (M = 0.47, SE = 0.16, p = 0.01). For condition 2, the constrast score (difference in means) between 1-week post and 14-week post was significantly different (M = 0.62, SE = 0.16, p = 0.00). So, this is basically another way of showing the main effect of time.

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_14wkpostEnd     -0.23931613
## 3                       Condition     -0.06821496
## 4 onewkpost_14wkpostEnd: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")

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

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”)
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.
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               Base         0 323        13.57223 3.910311 0.2175755 0.4280490
## 2               Base         1 347        12.90287 4.359070 0.2340071 0.4602555
## 3                End         0 344        12.65017 4.463040 0.2406310 0.4732981
## 4                End         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("Base"="1-Week Post", "End"="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_14wkpostEnd 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 Pr(>|t|)
## (Intercept)                      13.6835     0.3994 110.0831  34.263   <2e-16
## onewkpost_14wkpostEnd            -0.7267     0.3102 104.0786  -2.343    0.021
## Condition                        -0.6728     0.5609 108.8821  -1.200    0.233
## onewkpost_14wkpostEnd:Condition  -0.4486     0.4374 102.3356  -1.026    0.307
##                                    
## (Intercept)                     ***
## onewkpost_14wkpostEnd           *  
## Condition                          
## onewkpost_14wkpostEnd:Condition    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) on_14E Condtn
## onwkpst_14E -0.097              
## Condition   -0.712  0.069       
## onwkp_14E: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 [End]	-0.73	-1.33 – -0.12	0.019
Condition	-0.67	-1.77 – 0.43	0.230
onewkpost_14wkpost [End] * Condition	-0.45	-1.31 – 0.41	0.305
Random Effects
σ²	9.74
τ₀₀ _PID	7.06
τ₁₁ _{PID.onewkpost_14wkpostEnd}	1.92
ρ₀₁ _PID	0.29
ICC	0.48
N _PID	111
Observations	1328
Marginal R² / Conditional R²	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_14wkpostEnd     -0.21926308
## 3                       Condition     -0.07726937
## 4 onewkpost_14wkpostEnd: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?

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               Base         0 321        13.65368 3.782837 0.2111374 0.4153927
## 2               Base         1 347        12.90287 4.359070 0.2340071 0.4602555
## 3                End         0 332        13.05921 3.963221 0.2175100 0.4278763
## 4                End         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("Base"="1-Week Post", "End"="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 == 'Base')

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

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.

# Between Group Contrasts from the interaction model

#install.packages("emmeans")

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

est_means_out #Gives us all four simple means / SEs

##  onewkpost_14wkpost Condition emmean    SE  df lower.CL upper.CL
##  Base                       0   13.7 0.399 110     12.9     14.5
##  End                        0   13.0 0.481 108     12.0     13.9
##  Base                       1   13.0 0.394 108     12.2     13.8
##  End                        1   11.8 0.480 110     10.9     12.8
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

contrast_means_btw_out <- est_means_out %>% 
  emmeans::contrast(by = "onewkpost_14wkpost", simple = "Condition")

contrast_means_btw_out #Gives contrast stats for simple main effect at pre and then at post, with significance values

## onewkpost_14wkpost = Base:
##  contrast estimate   SE  df t.ratio p.value
##  0 effect    0.336 0.28 109  1.199  0.2330 
##  1 effect   -0.336 0.28 109 -1.199  0.2330 
## 
## onewkpost_14wkpost = End:
##  contrast estimate   SE  df t.ratio p.value
##  0 effect    0.561 0.34 109  1.650  0.1019 
##  1 effect   -0.561 0.34 109 -1.650  0.1019 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: fdr method for 2 tests

Interpreting the simple means (outputs from est_means)

Group 1 1-week post (M = 13.7, SE = 0.40) & end of semester (M = 13.0, SE = 0.48) Group 2 1-week post (M = 13.0, SE = 0.39) & end of semester (M = 11.8, SE = 0.48)

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

The contrast score (difference in means) between group 1 and group 2 was not significantly different at either 1-week post (M = -0.34, SE = 0.28, p = 0.25) or 14-weeks post (M = -0.56, SE = 0.34, p = 0.10).

Remaining questions - confirm how is this method dealing with the nested structure.

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 [End]	-0.07	-0.28 – 0.13	0.484
Condition	0.00	-0.29 – 0.30	0.993
onewkpost_14wkpost [End] * Condition	-0.15	-0.44 – 0.14	0.310
Random Effects
σ²	0.89
τ₀₀ _PID	0.47
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.30
ρ₀₁ _PID	-0.42
ICC	0.34
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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_14wkpostEnd   -0.1288100614
## 3                       Condition    0.0006049417
## 4 onewkpost_14wkpostEnd:Condition   -0.0649484143

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               Base         0 323            1.956048 1.221739 0.06797938
## 2               Base         1 347            1.926877 1.099753 0.05903781
## 3                End         0 344            1.871247 1.180213 0.06363285
## 4                End         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("Base"="1-Week Post", "End"="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

# Between Group Contrasts from the interaction model

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

est_means_neg #Gives us all four simple means / SEs

##  onewkpost_14wkpost Condition emmean    SE  df lower.CL upper.CL
##  Base                       0   1.95 0.107 110     1.74     2.16
##  End                        0   1.87 0.105 106     1.67     2.08
##  Base                       1   1.95 0.106 107     1.74     2.16
##  End                        1   1.72 0.106 110     1.51     1.93
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

contrast_means_btw_neg <- est_means_neg %>% 
  emmeans::contrast(by = "onewkpost_14wkpost", simple = "Condition")

contrast_means_btw_neg #Gives contrast stats for simple main effect at pre and then at post, with significance values

## onewkpost_14wkpost = Base:
##  contrast  estimate     SE  df t.ratio p.value
##  0 effect -0.000703 0.0752 108 -0.009  0.9926 
##  1 effect  0.000703 0.0752 108  0.009  0.9926 
## 
## onewkpost_14wkpost = End:
##  contrast  estimate     SE  df t.ratio p.value
##  0 effect  0.074723 0.0746 108  1.002  0.3186 
##  1 effect -0.074723 0.0746 108 -1.002  0.3186 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: fdr method for 2 tests

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 [End]	-0.33	-0.52 – -0.14	0.001
Condition	-0.09	-0.40 – 0.22	0.570
onewkpost_14wkpost [End] * Condition	-0.05	-0.31 – 0.22	0.722
Random Effects
σ²	0.65
τ₀₀ _PID	0.59
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.27
ρ₀₁ _PID	0.02
ICC	0.53
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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_14wkpostEnd     -0.29729505
## 3                       Condition     -0.03802438
## 4 onewkpost_14wkpostEnd: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 [End]	-0.25	-0.42 – -0.07	0.005
Condition	-0.14	-0.40 – 0.11	0.269
onewkpost_14wkpost [End] * Condition	-0.13	-0.38 – 0.12	0.301
Random Effects
σ²	0.60
τ₀₀ _PID	0.37
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.23
ρ₀₁ _PID	-0.11
ICC	0.43
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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_14wkpostEnd     -0.30462175
## 3                       Condition     -0.06944586
## 4 onewkpost_14wkpostEnd: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
##  Base                       0   3.22 0.0925 110     3.04     3.40
##  End                        0   2.97 0.1076 107     2.76     3.18
##  Base                       1   3.07 0.0912 107     2.89     3.26
##  End                        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 = Base:
##  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 = End:
##  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 [End]	-0.23	-0.38 – -0.09	0.002
Condition	-0.31	-0.54 – -0.08	0.009
onewkpost_14wkpost [End] * Condition	-0.11	-0.31 – 0.10	0.320
Random Effects
σ²	0.51
τ₀₀ _PID	0.30
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.13
ρ₀₁ _PID	-0.07
ICC	0.41
N _PID	110
Observations	1290
Marginal R² / Conditional R²	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_14wkpostEnd 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 Pr(>|t|)
## (Intercept)                       2.05936    0.15678 113.45533  13.136   <2e-16
## onewkpost_14wkpostEnd            -0.23918    0.16792 104.66484  -1.424   0.1573
## Condition                        -0.40727    0.21975 111.25272  -1.853   0.0665
## onewkpost_14wkpostEnd:Condition   0.03363    0.23742 104.35512   0.142   0.8876
##                                    
## (Intercept)                     ***
## onewkpost_14wkpostEnd              
## Condition                       .  
## onewkpost_14wkpostEnd:Condition    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) on_14E Condtn
## onwkpst_14E -0.305              
## Condition   -0.713  0.217       
## onwkp_14E: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 [End]	-0.24	-0.57 – 0.09	0.154
Condition	-0.41	-0.84 – 0.02	0.064
onewkpost_14wkpost [End] * Condition	0.03	-0.43 – 0.50	0.887
Random Effects
σ²	2.83
τ₀₀ _PID	0.86
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.60
ρ₀₁ _PID	0.07
ICC	0.30
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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_14wkpostEnd    -0.110367333
## 3                       Condition    -0.101168294
## 4 onewkpost_14wkpostEnd: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 [End]	-0.12	-0.41 – 0.18	0.440
Condition	-0.22	-0.55 – 0.11	0.182
onewkpost_14wkpost [End] * Condition	-0.05	-0.46 – 0.37	0.833
Random Effects
σ²	1.90
τ₀₀ _PID	0.44
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.55
ρ₀₁ _PID	0.04
ICC	0.28
N _PID	111
Observations	1253
Marginal R² / Conditional R²	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_14wkpostEnd     -0.08585512
## 3                       Condition     -0.06974864
## 4 onewkpost_14wkpostEnd: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 == 'Base')

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

###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 [End]	-0.14	-0.31 – 0.03	0.112
Condition	-0.02	-0.36 – 0.32	0.901
onewkpost_14wkpost [End] * Condition	-0.06	-0.31 – 0.18	0.600
Random Effects
σ²	0.64
τ₀₀ _PID	0.71
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.20
ρ₀₁ _PID	0.18
ICC	0.58
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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_14wkpostEnd    -0.138550060
## 3                       Condition    -0.008669619
## 4 onewkpost_14wkpostEnd: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 [End]	0.03	-0.06 – 0.11	0.535
Condition	0.04	-0.06 – 0.14	0.438
onewkpost_14wkpost [End] * Condition	-0.10	-0.22 – 0.02	0.097
Random Effects
σ²	0.22
τ₀₀ _PID	0.03
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.03
ρ₀₁ _PID	0.33
ICC	0.21
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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_14wkpostEnd     -0.04586647
## 3                       Condition      0.03740360
## 4 onewkpost_14wkpostEnd: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 [End]	0.03	-0.07 – 0.13	0.601
Condition	0.05	-0.02 – 0.12	0.188
onewkpost_14wkpost [End] * Condition	-0.09	-0.23 – 0.05	0.222
Random Effects
σ²	0.22
τ₀₀ _PID	0.00
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.07
ρ₀₁ _PID
N _PID	111
Observations	1253
Marginal R² / Conditional R²	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_14wkpostEnd     -0.03462212
## 3                       Condition      0.04880545
## 4 onewkpost_14wkpostEnd: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 [End]	-0.08	-0.24 – 0.08	0.320
Condition	-0.06	-0.27 – 0.15	0.585
onewkpost_14wkpost [End] * Condition	-0.07	-0.29 – 0.15	0.556
Random Effects
σ²	0.47
τ₀₀ _PID	0.25
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.19
ρ₀₁ _PID	-0.43
ICC	0.35
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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_14wkpostEnd     -0.13365482
## 3                       Condition     -0.03517613
## 4 onewkpost_14wkpostEnd: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
R² / R² 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
R² / R² 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 [End]	0.01	-0.07 – 0.08	0.865
Condition	0.06	-0.04 – 0.16	0.253
onewkpost_14wkpost [End] * Condition	-0.08	-0.19 – 0.02	0.126
Random Effects
σ²	0.22
τ₀₀ _PID	0.04
τ₁₁ _{PID.onewkpost_14wkpostEnd}	0.01
ρ₀₁ _PID	-0.17
ICC	0.16
N _PID	111
Observations	1345
Marginal R² / Conditional R²	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_14wkpostEnd     -0.06951643
## 3                       Condition      0.05901677
## 4 onewkpost_14wkpostEnd: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 [End]	-0.73	-1.33 – -0.12	0.019
Condition	-0.67	-1.77 – 0.43	0.230
onewkpost_14wkpost [End] * Condition	-0.45	-1.31 – 0.41	0.305
Random Effects
σ²	9.74
τ₀₀ _PID	7.06
τ₁₁ _{PID.onewkpost_14wkpostEnd}	1.92
ρ₀₁ _PID	0.29
ICC	0.48
N _PID	111
Observations	1328
Marginal R² / Conditional R²	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 == 'Base')

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

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
R² / R² 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
R² / R² 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
R² / R² adjusted	0.022 / 0.021

EMA Social Connection Scale Analysis

Megan Lipsett

Spring 2020

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

Correlation Matrix - All items

Items / Composites that will be analyzed

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

Negative Interactions

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

Re-shaping variables to factors

Interpreting the simple means (outputs from est_means)

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

Remaining questions - confirm how is this method dealing with the nested structure.

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

Effect Size

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

Plotting EMA_comp scores

Removing Outliers

Removing Outliers

Interpreting the simple means (outputs from est_means)

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

Remaining questions - confirm how is this method dealing with the nested structure.

Describing the data

MODELS FOR INDIVIDUAL ITEMS

Running models with outliers removed

Detecting Outliers

Removing outliers

Removing Outliers

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

Removing Outliers

Removing outliers

Removing Outliers

Removing outliers

Removing Outliers

Removing Outliers

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

EMA Social Connection Scale Analysis

Megan Lipsett

Spring 2020

SOCIAL CONNECTION SCALE ’

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

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

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.

Step 4: Impute series means for any remaining missing data

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.

Correlation Matrix - All items

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.

Write up- Internal validitiy of social connection

Checking for internal validity of positive social interaction composite.

Write up- Internal validity of positive social connection

Checking for internal validity of negative social interaction composite.

Write up- Internal validity of negative social connection

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.

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

Re-shaping variables to factors

Mixed models for EMA Social Connection Scales

POSITIVE, CLOSE, and CONNECTED SOCIAL CONNECTION COMPOSITE

Describing the data

Mixed model for Social Connection Composite (POST)

WRITE UP - Mixed Model Post Social Composite

Interpreting the simple means (outputs from est_means)

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

Remaining questions - confirm how is this method dealing with the nested structure.

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

Effect Size

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

Plotting EMA_comp scores

Removing Outliers

Removing Outliers

Mixed model for EMA Social Connection (outliers removed)

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

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

Plotting EMA_comp scores w/ outliers removed

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

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

Interpreting the simple means (outputs from est_means)

Interpreting the outputs from contrast_means_btw (Note: std error is the std err around the diff in means)

Remaining questions - confirm how is this method dealing with the nested structure.

///NEGATIVE SOCIAL INTERACTIONS COMPOSITE///

Describing the data

Mixed model for Negative Social Interactions Composite (POST)

Plotting EMA_neg_social_comp scores

MODELS FOR INDIVIDUAL ITEMS

NUMBER OF SOCIAL INTERACTIONS

Describing the data

Mixed model for Number of social interactions (POST)

NUMBER OF UNIQUE SOCIAL INTERACTIONS

Describing the data

Mixed model for Number of unique social interactions (POST)

Running models with outliers removed

Detecting Outliers

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

/// SOCIAL CLOSENESS / CONNECTION ///

Mixed model for social closeness (POST)

Removing outliers

Removing Outliers

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