andrea compassion and sleep
2023-08-29
Lib/function
library(psych)
library(nlme)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✔ ggplot2 3.4.2 ✔ purrr 1.0.1
## ✔ tibble 3.2.1 ✔ dplyr 1.1.2
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.4 ✔ forcats 1.0.0## Warning: package 'ggplot2' was built under R version 4.3.0## Warning: package 'tibble' was built under R version 4.3.1## Warning: package 'tidyr' was built under R version 4.2.2## Warning: package 'readr' was built under R version 4.2.2## Warning: package 'purrr' was built under R version 4.2.2## Warning: package 'dplyr' was built under R version 4.3.0## Warning: package 'stringr' was built under R version 4.3.0## Warning: package 'forcats' was built under R version 4.2.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ ggplot2::%+%() masks psych::%+%()
## ✖ ggplot2::alpha() masks psych::alpha()
## ✖ dplyr::collapse() masks nlme::collapse()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(ggplot2)
library(readit)
library(MuMIn)
ComptoSleep <- read.csv("~/andrea sleep and sc/data/compassionproceedingsleep.csv") #last self compassion before sleep
averageComptoSleep<-read.csv("~/andrea sleep and sc/data/AveragecompassionproceedingsleepA.csv") ##average for day
process_data <- function(Cortable) {
temp1 <- round(Cortable$rwg, 2)
temp1 <- temp1[order(rownames(temp1)), order(colnames(temp1))]
temp2 <- round(Cortable$pwg, 4)
temp2 <- temp2[order(rownames(temp2)), order(colnames(temp2))]
temp3 <- round(Cortable$rbg, 2)
temp3 <- temp3[order(rownames(temp3)), order(colnames(temp3))]
temp4 <- round(Cortable$pbg, 4)
temp4 <- temp4[order(rownames(temp4)), order(colnames(temp4))]
return(list(rwg = temp1, pwg = temp2, rbg = temp3, pbg = temp4))
}
information_theoretic_tables <- function(model) {
# Run the dredge function to perform multi-model selection
dr1 <- dredge(model)
# Get the model selection table
model_selection <- model.sel(dr1)
# Calculate the summed Akaike weights for each variable
var_importance <- sw(model.avg(dr1))
# Calculate the average parameter estimate across all models
avg_mod <- model.avg(dr1)
# Get the summary of the average model
avg_mod_summary <- summary(avg_mod)
# Get coefficients
coef <- coef(avg_mod, full = TRUE)
# Calculate confidence intervals for the average model
cis <- as.data.frame(confint(avg_mod, full=TRUE))
# Convert coefficients into a data frame and set row names
coef <- as.data.frame(coef)
rownames(coef) <- rownames(cis)
# Combine the coefficients with the confidence intervals
result <- cbind(coef, cis)
# Apply the round function to all numeric columns
result[] <- lapply(result, function(x) if(is.numeric(x)) round(x, 2) else x)
# Convert rownames into a column for variable names
result$vars <- rownames(result)
# Return a list containing the variable importance, model selection, result, and summary of the average model
return(list(var_importance = var_importance, model_selection = model_selection, result = result, avg_mod_summary = avg_mod_summary))
}R last
This is the time just before they went to sleep. i think we won’t focus on this but on the average of the day before they went to sleep
wg = within group correlation pwg = p value for that correlation
rbg=correlation between groups prb =p value of that correlation between groups
# Specify the variables you want to correlate
CorVariables <- c("SleepHours", "SleepQuality", "SleepRecovery", "last_SelfCompassion", "last_OtherCompassion","last_CurrentMood", "ID")
# Subset the data to include only the specified variables
subset_data <- ComptoSleep[, CorVariables]
# Calculate the correlations for each unique ID
Cortable <- statsBy(subset_data, "ID", cors = FALSE, cor = "cor", method = "pearson", use = "pairwise",
poly = FALSE, na.rm = TRUE, alpha = .05, minlength = 5, weights = NULL)
result <- process_data(Cortable)
result## $rwg
## last_CurrentMood.wg last_OtherCompassion.wg
## last_CurrentMood.wg 1.00 0.39
## last_OtherCompassion.wg 0.39 1.00
## last_SelfCompassion.wg 0.38 0.37
## SleepHours.wg 0.03 0.04
## SleepQuality.wg 0.06 0.08
## SleepRecovery.wg 0.16 0.13
## last_SelfCompassion.wg SleepHours.wg SleepQuality.wg
## last_CurrentMood.wg 0.38 0.03 0.06
## last_OtherCompassion.wg 0.37 0.04 0.08
## last_SelfCompassion.wg 1.00 0.01 0.05
## SleepHours.wg 0.01 1.00 0.45
## SleepQuality.wg 0.05 0.45 1.00
## SleepRecovery.wg 0.05 0.34 0.57
## SleepRecovery.wg
## last_CurrentMood.wg 0.16
## last_OtherCompassion.wg 0.13
## last_SelfCompassion.wg 0.05
## SleepHours.wg 0.34
## SleepQuality.wg 0.57
## SleepRecovery.wg 1.00
##
## $pwg
## last_CurrentMood.wg last_OtherCompassion.wg
## last_CurrentMood.wg 0.0000 0.0000
## last_OtherCompassion.wg 0.0000 0.0000
## last_SelfCompassion.wg 0.0000 0.0000
## SleepHours.wg 0.4611 0.2151
## SleepQuality.wg 0.0932 0.0207
## SleepRecovery.wg 0.0000 0.0001
## last_SelfCompassion.wg SleepHours.wg SleepQuality.wg
## last_CurrentMood.wg 0.0000 0.4611 0.0932
## last_OtherCompassion.wg 0.0000 0.2151 0.0207
## last_SelfCompassion.wg 0.0000 0.7069 0.1755
## SleepHours.wg 0.7069 0.0000 0.0000
## SleepQuality.wg 0.1755 0.0000 0.0000
## SleepRecovery.wg 0.1167 0.0000 0.0000
## SleepRecovery.wg
## last_CurrentMood.wg 0.0000
## last_OtherCompassion.wg 0.0001
## last_SelfCompassion.wg 0.1167
## SleepHours.wg 0.0000
## SleepQuality.wg 0.0000
## SleepRecovery.wg 0.0000
##
## $rbg
## last_CurrentMood.bg last_OtherCompassion.bg
## last_CurrentMood.bg 1.00 0.44
## last_OtherCompassion.bg 0.44 1.00
## last_SelfCompassion.bg 0.67 0.33
## SleepHours.bg 0.20 0.12
## SleepQuality.bg 0.53 0.23
## SleepRecovery.bg 0.56 0.27
## last_SelfCompassion.bg SleepHours.bg SleepQuality.bg
## last_CurrentMood.bg 0.67 0.20 0.53
## last_OtherCompassion.bg 0.33 0.12 0.23
## last_SelfCompassion.bg 1.00 0.16 0.48
## SleepHours.bg 0.16 1.00 0.45
## SleepQuality.bg 0.48 0.45 1.00
## SleepRecovery.bg 0.51 0.30 0.64
## SleepRecovery.bg
## last_CurrentMood.bg 0.56
## last_OtherCompassion.bg 0.27
## last_SelfCompassion.bg 0.51
## SleepHours.bg 0.30
## SleepQuality.bg 0.64
## SleepRecovery.bg 1.00
##
## $pbg
## last_CurrentMood.bg last_OtherCompassion.bg
## last_CurrentMood.bg 0.0000 0.0000
## last_OtherCompassion.bg 0.0000 0.0000
## last_SelfCompassion.bg 0.0000 0.0001
## SleepHours.bg 0.0217 0.1995
## SleepQuality.bg 0.0000 0.0100
## SleepRecovery.bg 0.0000 0.0027
## last_SelfCompassion.bg SleepHours.bg SleepQuality.bg
## last_CurrentMood.bg 0.0000 0.0217 0.00
## last_OtherCompassion.bg 0.0001 0.1995 0.01
## last_SelfCompassion.bg 0.0000 0.0749 0.00
## SleepHours.bg 0.0749 0.0000 0.00
## SleepQuality.bg 0.0000 0.0000 0.00
## SleepRecovery.bg 0.0000 0.0005 0.00
## SleepRecovery.bg
## last_CurrentMood.bg 0.0000
## last_OtherCompassion.bg 0.0027
## last_SelfCompassion.bg 0.0000
## SleepHours.bg 0.0005
## SleepQuality.bg 0.0000
## SleepRecovery.bg 0.0000R average
this is the average of the values on the day before sleep
CorVariables <- c("SleepHours","SleepQuality","SleepRecovery",
"SelfCompassion","OtherCompassion","CurrentMood",
"ID")
subset_data <- averageComptoSleep[, CorVariables]
# Calculate the correlations for each unique ID
Cortable <- statsBy(subset_data, "ID", cors = FALSE, cor = "cor", method = "pearson", use = "pairwise",
poly = FALSE, na.rm = TRUE, alpha = .05, minlength = 5, weights = NULL)
resultMean <- process_data(Cortable)
resultMean## $rwg
## CurrentMood.wg OtherCompassion.wg SelfCompassion.wg
## CurrentMood.wg 1.00 0.36 0.45
## OtherCompassion.wg 0.36 1.00 0.35
## SelfCompassion.wg 0.45 0.35 1.00
## SleepHours.wg 0.02 -0.04 -0.01
## SleepQuality.wg 0.01 0.05 0.07
## SleepRecovery.wg 0.08 0.10 0.11
## SleepHours.wg SleepQuality.wg SleepRecovery.wg
## CurrentMood.wg 0.02 0.01 0.08
## OtherCompassion.wg -0.04 0.05 0.10
## SelfCompassion.wg -0.01 0.07 0.11
## SleepHours.wg 1.00 0.45 0.34
## SleepQuality.wg 0.45 1.00 0.57
## SleepRecovery.wg 0.34 0.57 1.00
##
## $pwg
## CurrentMood.wg OtherCompassion.wg SelfCompassion.wg
## CurrentMood.wg 0.0000 0.0000 0.0000
## OtherCompassion.wg 0.0000 0.0000 0.0000
## SelfCompassion.wg 0.0000 0.0000 0.0000
## SleepHours.wg 0.6158 0.2897 0.7644
## SleepQuality.wg 0.6653 0.1269 0.0451
## SleepRecovery.wg 0.0236 0.0047 0.0019
## SleepHours.wg SleepQuality.wg SleepRecovery.wg
## CurrentMood.wg 0.6158 0.6654 0.0236
## OtherCompassion.wg 0.2897 0.1269 0.0047
## SelfCompassion.wg 0.7644 0.0451 0.0019
## SleepHours.wg 0.0000 0.0000 0.0000
## SleepQuality.wg 0.0000 0.0000 0.0000
## SleepRecovery.wg 0.0000 0.0000 0.0000
##
## $rbg
## CurrentMood.bg OtherCompassion.bg SelfCompassion.bg
## CurrentMood.bg 1.00 0.44 0.64
## OtherCompassion.bg 0.44 1.00 0.31
## SelfCompassion.bg 0.64 0.31 1.00
## SleepHours.bg 0.24 0.11 0.14
## SleepQuality.bg 0.58 0.23 0.47
## SleepRecovery.bg 0.62 0.25 0.52
## SleepHours.bg SleepQuality.bg SleepRecovery.bg
## CurrentMood.bg 0.24 0.58 0.62
## OtherCompassion.bg 0.11 0.23 0.25
## SelfCompassion.bg 0.14 0.47 0.52
## SleepHours.bg 1.00 0.45 0.30
## SleepQuality.bg 0.45 1.00 0.64
## SleepRecovery.bg 0.30 0.64 1.00
##
## $pbg
## CurrentMood.bg OtherCompassion.bg SelfCompassion.bg
## CurrentMood.bg 0.0000 0.0000 0.0000
## OtherCompassion.bg 0.0000 0.0000 0.0004
## SelfCompassion.bg 0.0000 0.0004 0.0000
## SleepHours.bg 0.0063 0.2148 0.1154
## SleepQuality.bg 0.0000 0.0089 0.0000
## SleepRecovery.bg 0.0000 0.0051 0.0000
## SleepHours.bg SleepQuality.bg SleepRecovery.bg
## CurrentMood.bg 0.0063 0.0000 0.0000
## OtherCompassion.bg 0.2148 0.0089 0.0051
## SelfCompassion.bg 0.1154 0.0000 0.0000
## SleepHours.bg 0.0000 0.0000 0.0005
## SleepQuality.bg 0.0000 0.0000 0.0000
## SleepRecovery.bg 0.0005 0.0000 0.0000Antecedent model
We are treating the compassion as an antecedent
random intercept
modelRIAVerage <- lme(scale(SleepRecovery) ~ scale(OtherCompassion) + scale(SelfCompassion)+ scale(CurrentMood)+
scale(SleepRecovery_lag),
random = ~ 1 | ID,
data = averageComptoSleep,
control = lmeControl(opt = 'optim'),
na.action = na.exclude)
summary(modelRIAVerage)## Linear mixed-effects model fit by REML
## Data: averageComptoSleep
## AIC BIC logLik
## 1900.009 1933.225 -943.0044
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.5912849 0.6246174
##
## Fixed effects: scale(SleepRecovery) ~ scale(OtherCompassion) + scale(SelfCompassion) + scale(CurrentMood) + scale(SleepRecovery_lag)
## Value Std.Error DF t-value p-value
## (Intercept) -0.00209281 0.05687039 725 -0.036800 0.9707
## scale(OtherCompassion) 0.04102188 0.03869160 725 1.060227 0.2894
## scale(SelfCompassion) 0.17673085 0.04644951 725 3.804795 0.0002
## scale(CurrentMood) 0.08190940 0.04057044 725 2.018943 0.0439
## scale(SleepRecovery_lag) 0.11868414 0.02926782 725 4.055107 0.0001
## Correlation:
## (Intr) sc(OC) sc(SC) sc(CM)
## scale(OtherCompassion) -0.001
## scale(SelfCompassion) 0.001 -0.176
## scale(CurrentMood) -0.001 -0.263 -0.400
## scale(SleepRecovery_lag) 0.001 0.007 -0.099 -0.178
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.21332295 -0.52272603 0.04782764 0.55520595 2.69692646
##
## Number of Observations: 855
## Number of Groups: 126random slope. (self-compasson as random )
modelRSLPAVerage <- lme(
scale(SleepRecovery) ~ scale(OtherCompassion) + scale(SelfCompassion) + scale(CurrentMood) + scale(SleepRecovery_lag),
random = ~ 1 + scale(SelfCompassion) | ID, ##only differs here
data = averageComptoSleep,
control = lmeControl(opt = 'optim'),
na.action = na.exclude
)
summary(modelRSLPAVerage)## Linear mixed-effects model fit by REML
## Data: averageComptoSleep
## AIC BIC logLik
## 1896.362 1939.069 -939.1809
##
## Random effects:
## Formula: ~1 + scale(SelfCompassion) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.5545538 (Intr)
## scale(SelfCompassion) 0.2491222 0.232
## Residual 0.6132156
##
## Fixed effects: scale(SleepRecovery) ~ scale(OtherCompassion) + scale(SelfCompassion) + scale(CurrentMood) + scale(SleepRecovery_lag)
## Value Std.Error DF t-value p-value
## (Intercept) -0.00320036 0.05646346 725 -0.056680 0.9548
## scale(OtherCompassion) 0.03695438 0.03857945 725 0.957877 0.3384
## scale(SelfCompassion) 0.18764323 0.05179737 725 3.622640 0.0003
## scale(CurrentMood) 0.08697639 0.04052341 725 2.146325 0.0322
## scale(SleepRecovery_lag) 0.11713432 0.02910284 725 4.024842 0.0001
## Correlation:
## (Intr) sc(OC) sc(SC) sc(CM)
## scale(OtherCompassion) 0.029
## scale(SelfCompassion) 0.074 -0.154
## scale(CurrentMood) -0.001 -0.255 -0.356
## scale(SleepRecovery_lag) 0.006 0.016 -0.089 -0.180
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.02684808 -0.51428261 0.05271255 0.55408636 2.75319273
##
## Number of Observations: 855
## Number of Groups: 126anova(modelRIAVerage,modelRSLPAVerage)## Model df AIC BIC logLik Test L.Ratio p-value
## modelRIAVerage 1 7 1900.009 1933.225 -943.0044
## modelRSLPAVerage 2 9 1896.362 1939.069 -939.1809 1 vs 2 7.646997 0.0219random slope (self-compassion and mood as random)
we can’t have too many random effects without model breaking, so we focus on significant IVS of interest
Adding mood doesn’t help so we just focus on self compassion as random, meaning it varies from person to person
modelRSLPAVerage2 <- lme(
scale(SleepRecovery) ~ scale(OtherCompassion) + scale(SelfCompassion) + scale(CurrentMood) + scale(SleepRecovery_lag),
random = ~ 1 + scale(SelfCompassion)+scale(CurrentMood) | ID, ##only differs here
data = averageComptoSleep,
control = lmeControl(opt = 'optim'),
na.action = na.exclude
)
summary(modelRSLPAVerage2)## Linear mixed-effects model fit by REML
## Data: averageComptoSleep
## AIC BIC logLik
## 1899.905 1956.848 -937.9524
##
## Random effects:
## Formula: ~1 + scale(SelfCompassion) + scale(CurrentMood) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.5543252 (Intr) sc(SC)
## scale(SelfCompassion) 0.2948517 0.277
## scale(CurrentMood) 0.1306583 -0.310 -0.678
## Residual 0.6081912
##
## Fixed effects: scale(SleepRecovery) ~ scale(OtherCompassion) + scale(SelfCompassion) + scale(CurrentMood) + scale(SleepRecovery_lag)
## Value Std.Error DF t-value p-value
## (Intercept) 0.01011670 0.05673791 725 0.178306 0.8585
## scale(OtherCompassion) 0.03630165 0.03893187 725 0.932440 0.3514
## scale(SelfCompassion) 0.19134902 0.05398260 725 3.544642 0.0004
## scale(CurrentMood) 0.08197907 0.04264113 725 1.922535 0.0549
## scale(SleepRecovery_lag) 0.11027309 0.02913853 725 3.784442 0.0002
## Correlation:
## (Intr) sc(OC) sc(SC) sc(CM)
## scale(OtherCompassion) 0.034
## scale(SelfCompassion) 0.103 -0.154
## scale(CurrentMood) -0.083 -0.236 -0.431
## scale(SleepRecovery_lag) 0.003 0.017 -0.089 -0.164
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.99828205 -0.51524756 0.04681509 0.54260603 2.75981230
##
## Number of Observations: 855
## Number of Groups: 126anova(modelRIAVerage,modelRSLPAVerage2)## Model df AIC BIC logLik Test L.Ratio p-value
## modelRIAVerage 1 7 1900.009 1933.225 -943.0044
## modelRSLPAVerage2 2 12 1899.905 1956.848 -937.9524 1 vs 2 10.10391 0.0723plot effects
fixed_effects <- fixef(modelRSLPAVerage)
# Extract random effects
random_effects <- ranef(modelRSLPAVerage)
# Convert the random effects to a data frame
random_effects_df <- as.data.frame(random_effects)
# Add a column with the random slopes for each person
random_effects_df$random_slope <- random_effects_df$`scale(SelfCompassion)`
# Add the fixed effect to each random slope to get the slope for each person
random_effects_df$person_slope <- fixed_effects[3] + random_effects_df$random_slope
# Calculate the mean of person_slope
mean_person_slope <- mean(random_effects_df$person_slope, na.rm = TRUE)
# Create the histogram
histogram <- ggplot(random_effects_df, aes(x=person_slope)) +
geom_histogram(color="black", fill="lightblue", bins=30) + # Use bins=30 as a starting point
geom_vline(aes(xintercept=mean_person_slope), color="red", linetype="dashed", linewidth=1) + # Replace size with linewidth
annotate("text", x = mean_person_slope, y= -Inf, label = paste("Mean =", round(mean_person_slope, 2)), vjust = -1.5) +
theme_minimal() +
labs(x="Person Slope", y="Count", title="Within person association: Self-compassion with sleep")
# Print the histogram
print(histogram)
Consequence model
We are treating compassion as a consequence of sleep
random intercept
we need two models, one for self compassion, one for other compassion
modelRIntAVerageSleepSelf <- lme(
scale(SelfCompassion) ~ scale(SleepHours_lag)+scale(SleepQuality_lag)+scale(SleepRecovery_lag)+scale(Mood_lag),
random = ~ 1 | ID,
data = averageComptoSleep,
control = lmeControl(opt = 'optim'),
na.action = na.exclude
)
summary(modelRIntAVerageSleepSelf)## Linear mixed-effects model fit by REML
## Data: averageComptoSleep
## AIC BIC logLik
## 1481.508 1514.6 -733.754
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.8486075 0.4502583
##
## Fixed effects: scale(SelfCompassion) ~ scale(SleepHours_lag) + scale(SleepQuality_lag) + scale(SleepRecovery_lag) + scale(Mood_lag)
## Value Std.Error DF t-value p-value
## (Intercept) -0.00015217 0.07719935 710 -0.001971 0.9984
## scale(SleepHours_lag) -0.04440681 0.02379206 710 -1.866455 0.0624
## scale(SleepQuality_lag) 0.03247852 0.02740041 710 1.185330 0.2363
## scale(SleepRecovery_lag) 0.10154829 0.02727870 710 3.722622 0.0002
## scale(Mood_lag) 0.01158144 0.02314234 710 0.500444 0.6169
## Correlation:
## (Intr) s(SH_) s(SQ_) s(SR_)
## scale(SleepHours_lag) 0.000
## scale(SleepQuality_lag) -0.001 -0.331
## scale(SleepRecovery_lag) 0.001 -0.105 -0.476
## scale(Mood_lag) -0.001 0.045 -0.062 -0.230
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.74837650 -0.51938058 0.08210379 0.50100857 3.77318097
##
## Number of Observations: 840
## Number of Groups: 126modelRIntAVerageSleepOther <- lme(
scale(OtherCompassion) ~ scale(SleepHours_lag)+scale(SleepQuality_lag)+scale(SleepRecovery_lag)+scale(Mood_lag),
random = ~ 1 | ID,
data = averageComptoSleep,
control = lmeControl(opt = 'optim'),
na.action = na.exclude
)
summary(modelRIntAVerageSleepOther)## Linear mixed-effects model fit by REML
## Data: averageComptoSleep
## AIC BIC logLik
## 1778.878 1811.97 -882.4392
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.8237174 0.5567971
##
## Fixed effects: scale(OtherCompassion) ~ scale(SleepHours_lag) + scale(SleepQuality_lag) + scale(SleepRecovery_lag) + scale(Mood_lag)
## Value Std.Error DF t-value p-value
## (Intercept) 0.00590783 0.07588511 710 0.0778523 0.9380
## scale(SleepHours_lag) -0.03274188 0.02912685 710 -1.1241132 0.2613
## scale(SleepQuality_lag) 0.02903727 0.03359692 710 0.8642839 0.3877
## scale(SleepRecovery_lag) 0.07118848 0.03340154 710 2.1312937 0.0334
## scale(Mood_lag) -0.01917601 0.02834120 710 -0.6766125 0.4989
## Correlation:
## (Intr) s(SH_) s(SQ_) s(SR_)
## scale(SleepHours_lag) 0.000
## scale(SleepQuality_lag) -0.001 -0.333
## scale(SleepRecovery_lag) 0.001 -0.103 -0.476
## scale(Mood_lag) -0.001 0.044 -0.066 -0.233
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.63525351 -0.37609714 0.09218939 0.46440134 3.18502816
##
## Number of Observations: 840
## Number of Groups: 126random slopes
sleep recovery looks like the big winner, so we make that random
modelRSLPAVerageSleepSelf <- lme(
scale(SelfCompassion) ~ scale(SleepHours_lag)+scale(SleepQuality_lag)+scale(SleepRecovery_lag)+scale(Mood_lag),
random = ~ 1+scale(SleepRecovery_lag) | ID,
data = averageComptoSleep,
control = lmeControl(opt = 'optim'),
na.action = na.exclude
)
summary(modelRSLPAVerageSleepSelf )## Linear mixed-effects model fit by REML
## Data: averageComptoSleep
## AIC BIC logLik
## 1481.168 1523.715 -731.5842
##
## Random effects:
## Formula: ~1 + scale(SleepRecovery_lag) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.8453021 (Intr)
## scale(SleepRecovery_lag) 0.1181127 0.035
## Residual 0.4412491
##
## Fixed effects: scale(SelfCompassion) ~ scale(SleepHours_lag) + scale(SleepQuality_lag) + scale(SleepRecovery_lag) + scale(Mood_lag)
## Value Std.Error DF t-value p-value
## (Intercept) -0.00379010 0.07714758 710 -0.049128 0.9608
## scale(SleepHours_lag) -0.04510621 0.02434129 710 -1.853074 0.0643
## scale(SleepQuality_lag) 0.03577819 0.02768292 710 1.292428 0.1966
## scale(SleepRecovery_lag) 0.10594844 0.03004930 710 3.525821 0.0004
## scale(Mood_lag) 0.01606801 0.02321520 710 0.692133 0.4891
## Correlation:
## (Intr) s(SH_) s(SQ_) s(SR_)
## scale(SleepHours_lag) -0.002
## scale(SleepQuality_lag) 0.000 -0.327
## scale(SleepRecovery_lag) 0.009 -0.086 -0.444
## scale(Mood_lag) 0.001 0.042 -0.061 -0.207
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.84404227 -0.52629427 0.07073492 0.51701144 3.83857487
##
## Number of Observations: 840
## Number of Groups: 126modelRSLPAVerageSleepOther <- lme(
scale(OtherCompassion) ~ scale(SleepHours_lag)+scale(SleepQuality_lag)+scale(SleepRecovery_lag)+scale(Mood_lag),
random = ~ 1 +scale(SleepRecovery_lag) | ID,
data = averageComptoSleep,
control = lmeControl(opt = 'optim'),
na.action = na.exclude
)
summary(modelRSLPAVerageSleepOther)## Linear mixed-effects model fit by REML
## Data: averageComptoSleep
## AIC BIC logLik
## 1766.827 1809.374 -874.4134
##
## Random effects:
## Formula: ~1 + scale(SleepRecovery_lag) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.8046866 (Intr)
## scale(SleepRecovery_lag) 0.1798071 -0.373
## Residual 0.5401469
##
## Fixed effects: scale(OtherCompassion) ~ scale(SleepHours_lag) + scale(SleepQuality_lag) + scale(SleepRecovery_lag) + scale(Mood_lag)
## Value Std.Error DF t-value p-value
## (Intercept) 0.01869335 0.07469028 710 0.2502782 0.8024
## scale(SleepHours_lag) -0.03856606 0.02989479 710 -1.2900596 0.1974
## scale(SleepQuality_lag) 0.01064651 0.03373177 710 0.3156225 0.7524
## scale(SleepRecovery_lag) 0.08004566 0.03814525 710 2.0984440 0.0362
## scale(Mood_lag) -0.02193227 0.02834091 710 -0.7738732 0.4393
## Correlation:
## (Intr) s(SH_) s(SQ_) s(SR_)
## scale(SleepHours_lag) -0.004
## scale(SleepQuality_lag) 0.000 -0.322
## scale(SleepRecovery_lag) -0.164 -0.078 -0.431
## scale(Mood_lag) -0.001 0.042 -0.063 -0.202
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.1170056 -0.3471619 0.1095493 0.4539483 3.3489581
##
## Number of Observations: 840
## Number of Groups: 126test for differences between random intercept and slope
anova(modelRIntAVerageSleepSelf ,modelRSLPAVerageSleepSelf) ##no random effects## Model df AIC BIC logLik Test L.Ratio
## modelRIntAVerageSleepSelf 1 7 1481.508 1514.600 -733.7540
## modelRSLPAVerageSleepSelf 2 9 1481.168 1523.715 -731.5842 1 vs 2 4.339466
## p-value
## modelRIntAVerageSleepSelf
## modelRSLPAVerageSleepSelf 0.1142anova(modelRIntAVerageSleepOther,modelRSLPAVerageSleepOther) #3biggish random effect here## Model df AIC BIC logLik Test L.Ratio
## modelRIntAVerageSleepOther 1 7 1778.878 1811.970 -882.4392
## modelRSLPAVerageSleepOther 2 9 1766.827 1809.374 -874.4134 1 vs 2 16.05145
## p-value
## modelRIntAVerageSleepOther
## modelRSLPAVerageSleepOther 3e-04plot random effect of sleep on other compassion
fixed_effects <- fixef(modelRSLPAVerageSleepOther)
# Extract random effects
random_effects <- ranef(modelRSLPAVerageSleepOther)
# Convert the random effects to a data frame
random_effects_df <- as.data.frame(random_effects)
# Add a column with the random slopes for each person
random_effects_df$random_slope <- random_effects_df$`scale(SleepRecovery_lag)`
# Add the fixed effect to each random slope to get the slope for each person
random_effects_df$person_slope <- fixed_effects[4] + random_effects_df$random_slope
# Calculate the mean of person_slope
mean_person_slope <- mean(random_effects_df$person_slope, na.rm = TRUE)
# Create the histogram
histogram <- ggplot(random_effects_df, aes(x=person_slope)) +
geom_histogram(color="black", fill="lightblue", bins=30) + # Use bins=30 as a starting point
geom_vline(aes(xintercept=mean_person_slope), color="red", linetype="dashed", linewidth=1) + # Replace size with linewidth
annotate("text", x = mean_person_slope, y= -Inf, label = paste("Mean =", round(mean_person_slope, 2)), vjust = -1.5) +
theme_minimal() +
labs(x="Person Slope", y="Count", title="Within person association: Self-compassion with sleep")
# Print the histogram
print(histogram)