Danika parenting and nonattachment
Lib/function
###libraries and functions#####
library(psych)
library(readit)
library(nlme)
library(MuMIn)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.0##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.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
## ✔ purrr 1.0.1## 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(Hmisc)##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following object is masked from 'package:psych':
##
## describe## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(DT)## Warning: package 'DT' was built under R version 4.3.0rm(list=ls()) ##clean up previous files
FullDanikaData <- readit("~/danika parenting and compassion/data/ACSReducedLong.csv") ##long just to #check copy mean## File guessed to be CSV ("~/danika parenting and compassion/data/ACSReducedLong.csv")## Rows: 7041 Columns: 24
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): sdem.ethd.t03, GenderFem
## dbl (22): sdem.age.t03, sdem.age.t04, sdem.age.t05, sdem.mar.t03, sdem.marot...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.ImputedDanikaData<-readit("~/danika parenting and compassion/data/ACS_imputed.csv") #long## File guessed to be CSV ("~/danika parenting and compassion/data/ACS_imputed.csv")
## Rows: 7041 Columns: 8── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): GenderFem
## dbl (7): ID, Year, Pmon, NonAttach, Pauth, PControl, PSupport
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.acsClean <- read.csv("~/danika parenting and compassion/data/ACSReducedWide.csv") #wide
n_distinct(FullDanikaData$ID) ## [1] 2347n_distinct(ImputedDanikaData$ID)## [1] 2347n_distinct(acsClean$ID)## [1] 2347#info theoretic modelling; no iac cuttof###
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))
}
##information theoretic modeling: requires less than 7
# information_theoretic_tables <- function(model) {
# # Run dredge function for multi-model selection
# dr1 <- dredge(model)
#
# # Calculate the best (smallest) AIC value
# best_aic <- min(dr1$AIC)
#
# # Calculate ΔAIC for each model
# delta_aic <- dr1$AIC - best_aic
#
# # Filter models based on ΔAIC
# filtered_models <- dr1[delta_aic < 7, ]
#
# # Check for multiple models
# if (nrow(filtered_models) <= 1) {
# stop("Not enough models with ΔAIC < 7 for model averaging")
# }
#
# # Perform model selection on filtered models
# model_selection <- model.sel(filtered_models)
#
# # Calculate summed Akaike weights for each variable
# var_importance <- sw(model.avg(filtered_models))
#
# # Calculate average parameter estimate across filtered models
# avg_mod <- model.avg(filtered_models)
#
# # Get summary of average model
# avg_mod_summary <- summary(avg_mod)
#
# # Get coefficients
# coef <- coef(avg_mod, full = TRUE)
#
# # Calculate confidence intervals
# cis <- as.data.frame(confint(avg_mod, full=TRUE))
#
# # Convert coefficients to data frame
# coef <- as.data.frame(coef)
# rownames(coef) <- rownames(cis)
#
# # Combine coefficients and confidence intervals
# result <- cbind(coef, cis)
#
# # Round all numeric columns
# result[] <- lapply(result, function(x) if(is.numeric(x)) round(x, 2) else x)
#
# # Add variable names column
# result$vars <- rownames(result)
#
# # Return results
# return(list(var_importance = var_importance, model_selection = model_selection, result = result, avg_mod_summary = avg_mod_summary))
# }
options(scipen = 999) ##turn off scientific notation
# Robust control settings
ctrl <- lmeControl(opt="optim", msMaxIter = 500, niterEM=100, msMaxEval = 500, returnObject = TRUE, tolerance = 1e-06)
##function to calculate betweeen and within correlations####
calculateCorrelations <- function(df) {
# Calculate correlations
Cortable <- statsBy(df, "ID", cors = FALSE, cor = "cor", method = "pearson",
use = "pairwise", poly = FALSE, na.rm = TRUE,
alpha = 0.05, minlength = 5, weights = NULL)
# Round and order correlation matrices
temp1 <- round(Cortable$rwg, 2)
temp2 <- round(Cortable$pwg, 2)
temp3 <- round(Cortable$rbg, 2)
temp4 <- round(Cortable$pbg, 2)
# Create the list of items
result <- list(rwg = temp1, pwg = temp2, rbg = temp3, pbg = temp4)
return(result)
}
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))
}Missing
Danika, there are small, significant links between missingess and the variables. Generally, more missingness means less non attachment adn more negative parenting, but effects are small. acknowledge this in discussion section that your sample might have had skewed slightly because the kids with worse parenting experience were slighlty less likely to be present for testing…
# Calculate missingness by Year for each variable
variables_to_check <- c("NonAttach", "Pauth", "PControl", "Pmon", "PSupport")
table_results_list <- lapply(variables_to_check, function(var_name) {
FullDanikaData %>%
dplyr::group_by(Year) %>%
dplyr::summarize(
non_missing_values = sum(!is.na(.data[[var_name]])),
missing_values = sum(is.na(.data[[var_name]])),
percentage_non_missing = round((non_missing_values / (non_missing_values + missing_values)) * 100, 2)
) %>%
tibble::add_column(variable = var_name, .before = 1)
})
final_table_result <- do.call(rbind, table_results_list)
final_table_result## # A tibble: 15 × 5
## variable Year non_missing_values missing_values percentage_non_missing
## <chr> <dbl> <int> <int> <dbl>
## 1 NonAttach 10 1906 441 81.2
## 2 NonAttach 11 1664 683 70.9
## 3 NonAttach 12 1601 746 68.2
## 4 Pauth 10 1955 392 83.3
## 5 Pauth 11 1703 644 72.6
## 6 Pauth 12 1624 723 69.2
## 7 PControl 10 1953 394 83.2
## 8 PControl 11 1701 646 72.5
## 9 PControl 12 1628 719 69.4
## 10 Pmon 10 1952 395 83.2
## 11 Pmon 11 1704 643 72.6
## 12 Pmon 12 1630 717 69.4
## 13 PSupport 10 1954 393 83.3
## 14 PSupport 11 1695 652 72.2
## 15 PSupport 12 1628 719 69.4# Calculate missing points for each ID
summary_data <- data.frame()
for (var in variables_to_check) {
missing_var_name <- paste0("Missing_", var)
temp_data <- FullDanikaData %>%
group_by(ID) %>%
summarise(!!missing_var_name := sum(is.na(.data[[var]])))
if (nrow(summary_data) == 0) {
summary_data <- temp_data
} else {
summary_data <- left_join(summary_data, temp_data, by = "ID")
}
}
# Merge summary_data into FullDanikaData
FullDanikaData <- left_join(FullDanikaData, summary_data, by = "ID")
# Check for correlations between missing_vars and original_vars
missing_vars <- paste0("Missing_", variables_to_check)
original_vars <- variables_to_check
correlation_results <- data.frame()
for (i in 1:length(missing_vars)) {
# Check if columns exist
if (all(c(missing_vars[i], original_vars[i]) %in% names(FullDanikaData))) {
result <- rcorr(as.matrix(FullDanikaData[, c(missing_vars[i], original_vars[i])]), type = "pearson")
temp_df <- data.frame(
Variable1 = missing_vars[i],
Variable2 = original_vars[i],
Correlation = round(result$r[1, 2], 2),
Pvalue = round(result$P[1, 2], 2)
)
correlation_results <- rbind(correlation_results, temp_df)
} else {
message(paste("Skipped: Columns", missing_vars[i], "or", original_vars[i], "do not exist."))
}
}
# Display correlation results
print(correlation_results)## Variable1 Variable2 Correlation Pvalue
## 1 Missing_NonAttach NonAttach -0.06 0.00
## 2 Missing_Pauth Pauth -0.07 0.00
## 3 Missing_PControl PControl 0.05 0.00
## 4 Missing_Pmon Pmon -0.11 0.00
## 5 Missing_PSupport PSupport -0.04 0.01Descript
Age
mean(acsClean$sdem.age.t03, na.rm = TRUE)## [1] 15.64893sd(acsClean$sdem.age.t03, na.rm = TRUE)## [1] 0.4331287 mean(acsClean$sdem.age.t04, na.rm = TRUE)## [1] 16.58216sd(acsClean$sdem.age.t04, na.rm = TRUE)## [1] 0.4585875mean(acsClean$sdem.age.t05, na.rm = TRUE)## [1] 17.72658 sd(acsClean$sdem.age.t05, na.rm = TRUE)## [1] 0.3963813table(acsClean$GenderFem)##
## Female Male
## 1184 1157marital stats status 1) Married 2) Separated 3) Divorced 4)
other 1 = living together, 2 = mother deceased, 3 = father deceased, 4 = never together, 5 = unknown
table (acsClean$sdem.mar.t03)##
## 1 2 3 4
## 1433 177 248 80table(acsClean$sdem.marot.t03)##
## 1 2 3 4 5
## 46 9 11 3 7*********employment code this taken in previous year year 9 so lots of misisng 1 unemployed , 2 working, 3 intermediate, 4 salariate
not enough unemployed # SES. Child social class was based on parents’ # occupation recoded into the three classes (working, intermediate, and salariat classes) #using the Erickson–Goldthorpe–Portocarero (EGP) schema (Erikson, 1984). As per Erikson (1984) and Morgan, Spiller and Todd (2013), the ‘highest’ (most prestigious) # occupation of the child’s mother or father is used.
table_counts <- table(acsClean$`ses.totscore.t02`)
table_percentages <- prop.table(table_counts) * 100
result_table <- data.frame(
Count = table_counts,
Percentage = table_percentages
)
result_table## Count.Var1 Count.Freq Percentage.Var1 Percentage.Freq
## 1 1 8 1 0.5490734
## 2 2 110 2 7.5497598
## 3 3 438 3 30.0617708
## 4 4 901 4 61.8393960ethnicity 1 = aboroginal; 2 = torre strait islander; 3=caribean, 4=african; 5black american
6 caucasion 7 new zealand 8 american 9 european 10 chinese 11 japenese #12 vietnemes 17 other
table (acsClean$sdem.eth.t03) ##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 14 9 2 9 2 1036 13 6 155 11 5 8 4 18 3 4
## 17
## 645#table(acsClean$sdem.ethd.t03)Trends
means and standard devations. not much going on here at the mean level
# Calculate mean, standard deviation, and standard error
summary_table <- FullDanikaData %>%
group_by(Year) %>%
summarise(
NonAttach_Mean = mean(NonAttach, na.rm = TRUE),
NonAttach_SD = sd(NonAttach, na.rm = TRUE),
NonAttach_SE = sd(NonAttach, na.rm = TRUE) / sqrt(n()),
Pauth_Mean = mean(Pauth, na.rm = TRUE),
Pauth_SD = sd(Pauth, na.rm = TRUE),
Pauth_SE = sd(Pauth, na.rm = TRUE) / sqrt(n()),
PControl_Mean = mean(PControl, na.rm = TRUE),
PControl_SD = sd(PControl, na.rm = TRUE),
PControl_SE = sd(PControl, na.rm = TRUE) / sqrt(n()),
Pmon_Mean = mean(Pmon, na.rm = TRUE),
Pmon_SD = sd(Pmon, na.rm = TRUE),
Pmon_SE = sd(Pmon, na.rm = TRUE) / sqrt(n()),
PSupport_Mean = mean(PSupport, na.rm = TRUE),
PSupport_SD = sd(PSupport, na.rm = TRUE),
PSupport_SE = sd(PSupport, na.rm = TRUE) / sqrt(n())
)
# Display the summary table
datatable(summary_table)creating a nice ggplot table
# Reshape the data to a long format
long_data <- summary_table %>%
pivot_longer(
cols = -Year,
names_to = c("Variable", "Statistic"),
names_sep = "_",
values_to = "Value"
) %>%
filter(Statistic == "Mean")
# Create the plot
ggplot(long_data, aes(x = Year, y = Value, color = Variable)) +
geom_line() +
geom_point() +
labs(
title = "Mean Values by Year",
x = "Year",
y = "Mean Value"
) +
theme_minimal() +
scale_color_brewer(palette = "Set1")
Correl
CorImpute<- dplyr::select(ImputedDanikaData, -c(Year,GenderFem))
resultImpute<- calculateCorrelations(CorImpute)
process_data(resultImpute)## $rwg
## NonAttach.wg Pauth.wg PControl.wg Pmon.wg PSupport.wg
## NonAttach.wg 1.00 0.15 -0.10 0.11 0.19
## Pauth.wg 0.15 1.00 -0.20 0.24 0.37
## PControl.wg -0.10 -0.20 1.00 -0.20 -0.34
## Pmon.wg 0.11 0.24 -0.20 1.00 0.27
## PSupport.wg 0.19 0.37 -0.34 0.27 1.00
##
## $pwg
## NonAttach.wg Pauth.wg PControl.wg Pmon.wg PSupport.wg
## NonAttach.wg 0 0 0 0 0
## Pauth.wg 0 0 0 0 0
## PControl.wg 0 0 0 0 0
## Pmon.wg 0 0 0 0 0
## PSupport.wg 0 0 0 0 0
##
## $rbg
## NonAttach.bg Pauth.bg PControl.bg Pmon.bg PSupport.bg
## NonAttach.bg 1.00 0.36 -0.33 0.29 0.39
## Pauth.bg 0.36 1.00 -0.47 0.48 0.68
## PControl.bg -0.33 -0.47 1.00 -0.39 -0.61
## Pmon.bg 0.29 0.48 -0.39 1.00 0.49
## PSupport.bg 0.39 0.68 -0.61 0.49 1.00
##
## $pbg
## NonAttach.bg Pauth.bg PControl.bg Pmon.bg PSupport.bg
## NonAttach.bg 0 0 0 0 0
## Pauth.bg 0 0 0 0 0
## PControl.bg 0 0 0 0 0
## Pmon.bg 0 0 0 0 0
## PSupport.bg 0 0 0 0 0Contempor model
One variable at a time Nonattachment
we will only report the random slope as this is always the significant model
# Run random intercept model. what is the pooled effect/average effect
#of support on nonattaachment, within person (nested in id)
modelRI_Psupport <- lme(scale(NonAttach) ~ scale(PSupport), random = ~ 1 | ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
# Run random slope model.is the pooled effect the same for all people
modelRSslp_Psupport <- lme(scale(NonAttach) ~ scale(PSupport),
random = ~ 1 + scale(PSupport) | ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
anovaPsupport<-anova(modelRI_Psupport,modelRSslp_Psupport)
sumPsupport<-summary(modelRSslp_Psupport)
# Model using Pauth
modelRI_Pauth <- lme(scale(NonAttach) ~ scale(Pauth), random = ~ 1 | ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
modelRSslp_Pauth <- lme(scale(NonAttach) ~ scale(Pauth),
random = ~ 1 + scale(Pauth) | ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
anovaPauth <- anova(modelRI_Pauth, modelRSslp_Pauth)
sumPauth <- summary(modelRSslp_Pauth)
# Model using PControl
modelRI_PControl <- lme(scale(NonAttach) ~ scale(PControl), random = ~ 1 | ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
modelRSslp_PControl <- lme(scale(NonAttach) ~ scale(PControl),
random = ~ 1 + scale(PControl) | ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
anovaPControl <- anova(modelRI_PControl, modelRSslp_PControl)
sumPControl <- summary(modelRSslp_PControl)
# Model using Pmon
modelRI_Pmon <- lme(scale(NonAttach) ~ scale(Pmon), random = ~ 1 | ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
modelRSslp_Pmon <- lme(scale(NonAttach) ~ scale(Pmon),
random = ~ 1 + scale(Pmon) | ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
anovaPmon <- anova(modelRI_Pmon, modelRSslp_Pmon)
sumPmon <- summary(modelRSslp_Pmon)
anovaPauth## Model df AIC BIC logLik Test L.Ratio p-value
## modelRI_Pauth 1 5 16079.09 16113.39 -8034.545
## modelRSslp_Pauth 2 7 15774.93 15822.94 -7880.465 1 vs 2 308.1605 <.0001sumPauth## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 15774.93 15822.94 -7880.465
##
## Random effects:
## Formula: ~1 + scale(Pauth) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.6280905 (Intr)
## scale(Pauth) 0.4075250 0.066
## Residual 0.6176331
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.4056105
## Fixed effects: scale(NonAttach) ~ scale(Pauth)
## Value Std.Error DF t-value p-value
## (Intercept) 0.00253146 0.01713064 4693 0.147774 0.8825
## scale(Pauth) 0.24420181 0.01587488 4693 15.382912 0.0000
## Correlation:
## (Intr)
## scale(Pauth) -0.021
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.12876639 -0.45457310 -0.04139001 0.42143655 4.15854487
##
## Number of Observations: 7041
## Number of Groups: 2347anovaPControl## Model df AIC BIC logLik Test L.Ratio
## modelRI_PControl 1 5 16208.49 16242.78 -8099.243
## modelRSslp_PControl 2 7 15837.91 15885.92 -7911.954 1 vs 2 374.5776
## p-value
## modelRI_PControl
## modelRSslp_PControl <.0001sumPControl## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 15837.91 15885.92 -7911.954
##
## Random effects:
## Formula: ~1 + scale(PControl) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.6069852 (Intr)
## scale(PControl) 0.4430597 -0.162
## Residual 0.6332992
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.4474301
## Fixed effects: scale(NonAttach) ~ scale(PControl)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0063335 0.01730960 4693 0.365895 0.7145
## scale(PControl) -0.2029037 0.01596811 4693 -12.706806 0.0000
## Correlation:
## (Intr)
## scale(PControl) -0.025
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.45025783 -0.46330416 -0.04867964 0.42576568 4.16050111
##
## Number of Observations: 7041
## Number of Groups: 2347anovaPmon## Model df AIC BIC logLik Test L.Ratio p-value
## modelRI_Pmon 1 5 16223.30 16257.59 -8106.649
## modelRSslp_Pmon 2 7 15892.18 15940.20 -7939.092 1 vs 2 335.1154 <.0001sumPmon## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 15892.18 15940.2 -7939.092
##
## Random effects:
## Formula: ~1 + scale(Pmon) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.6258024 (Intr)
## scale(Pmon) 0.4148214 0.185
## Residual 0.6359717
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.4384988
## Fixed effects: scale(NonAttach) ~ scale(Pmon)
## Value Std.Error DF t-value p-value
## (Intercept) 0.01168255 0.01748687 4693 0.668075 0.5041
## scale(Pmon) 0.19935180 0.01593845 4693 12.507606 0.0000
## Correlation:
## (Intr)
## scale(Pmon) 0.011
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.47212747 -0.45635660 -0.04576391 0.43677493 3.39034026
##
## Number of Observations: 7041
## Number of Groups: 2347anovaPsupport## Model df AIC BIC logLik Test L.Ratio
## modelRI_Psupport 1 5 15922.17 15956.47 -7956.085
## modelRSslp_Psupport 2 7 15539.86 15587.87 -7762.929 1 vs 2 386.3133
## p-value
## modelRI_Psupport
## modelRSslp_Psupport <.0001sumPsupport## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 15539.86 15587.87 -7762.929
##
## Random effects:
## Formula: ~1 + scale(PSupport) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.5564232 (Intr)
## scale(PSupport) 0.4264815 0.228
## Residual 0.6451756
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.4751141
## Fixed effects: scale(NonAttach) ~ scale(PSupport)
## Value Std.Error DF t-value p-value
## (Intercept) -0.0027026 0.01670034 4693 -0.161829 0.8714
## scale(PSupport) 0.2969063 0.01583378 4693 18.751443 0.0000
## Correlation:
## (Intr)
## scale(PSupport) 0.024
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.86681882 -0.47366926 -0.03197217 0.45343986 3.85362471
##
## Number of Observations: 7041
## Number of Groups: 2347Plot the social support effects only, for sake of parsemony
##nicer plot
# Extract fixed effects
fixed_effects <- fixef(modelRSslp_Psupport)
# Extract random effects
random_effects <- ranef(modelRSslp_Psupport)
# 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(PSupport)`
# Add the fixed effect to each random slope to get the slope for each person
random_effects_df$person_slope <- fixed_effects[2] + 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", size=1) +
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: PSupport with Nonattachment")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# Print the histogram
print(histogram)
Full model with all variables. information theoretic modelling will run all possible models and determine importance of each variable. i can’t get it to run when estimating all random slopes together so we have to leave those out
support is the most important because confidence intervals don’t overlap
modelRIFull <- lme(scale(NonAttach) ~ scale(PSupport)+scale(Pmon)+scale(Pauth)+ scale(PControl),
random = ~ 1 | ID,
correlation = corAR1(),
control=ctrl,
data = ImputedDanikaData)
summary(modelRIFull)$tTable## Value Std.Error DF t-value
## (Intercept) 0.01178830 0.01655959 4690 0.7118716
## scale(PSupport) 0.17847320 0.01526939 4690 11.6882968
## scale(Pmon) 0.07323004 0.01302478 4690 5.6223618
## scale(Pauth) 0.10864715 0.01407776 4690 7.7176458
## scale(PControl) -0.06367897 0.01336532 4690 -4.7644919
## p-value
## (Intercept) 0.4765797123564615711543979159614536911
## scale(PSupport) 0.0000000000000000000000000000003938142
## scale(Pmon) 0.0000000199281941631764342921213055285
## scale(Pauth) 0.0000000000000143836731322315493636074
## scale(PControl) 0.0000019504831456704832820666598719450contemporanous_models<-information_theoretic_tables(modelRIFull)## Warning in dredge(model): comparing models fitted by REML## Fixed term is "(Intercept)"print(contemporanous_models$var_importance)## scale(PSupport) scale(Pauth) scale(Pmon) scale(PControl)
## Sum of weights: 1 1 1 1
## N containing models: 8 8 8 8print(contemporanous_models$model_selection)## Global model call: lme.formula(fixed = scale(NonAttach) ~ scale(PSupport) + scale(Pmon) +
## scale(Pauth) + scale(PControl), data = ImputedDanikaData,
## random = ~1 | ID, correlation = corAR1(), control = ctrl)
## ---
## Model selection table
## (Int) scl(Pth) scl(PCn) scl(Pmn) scl(PSp) df logLik AICc delta
## 16 0.01179 0.1086 -0.06368 0.07323 0.1785 8 -7891.944 15799.9 0.00
## 14 0.01107 0.1138 0.08002 0.2025 7 -7899.814 15813.6 13.74
## 12 0.01156 0.1230 -0.07194 0.1936 7 -7904.260 15822.5 22.63
## 10 0.01093 0.1304 0.2226 6 -7915.361 15842.7 42.83
## 15 0.01191 -0.07169 0.09145 0.2234 7 -7918.233 15850.5 50.57
## 13 0.01114 0.10010 0.2528 6 -7929.081 15870.2 70.27
## 11 0.01170 -0.08365 0.2504 6 -7939.923 15891.9 91.95
## 9 0.01099 0.2883 5 -7956.085 15922.2 122.27
## 8 0.01348 0.1714 -0.11560 0.10010 7 -7956.269 15926.6 126.65
## 4 0.01361 0.1991 -0.13340 6 -7982.614 15977.2 177.33
## 6 0.01287 0.1987 0.12130 6 -7993.510 15999.0 199.12
## 2 0.01300 0.2382 5 -8034.545 16079.1 279.19
## 7 0.01459 -0.15380 0.14590 6 -8036.675 16085.4 285.45
## 3 0.01514 -0.19060 5 -8099.243 16208.5 408.59
## 5 0.01407 0.18470 5 -8106.649 16223.3 423.40
## 1 0.01497 4 -8212.176 16432.4 632.45
## weight
## 16 0.999
## 14 0.001
## 12 0.000
## 10 0.000
## 15 0.000
## 13 0.000
## 11 0.000
## 9 0.000
## 8 0.000
## 4 0.000
## 6 0.000
## 2 0.000
## 7 0.000
## 3 0.000
## 5 0.000
## 1 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | IDprint(contemporanous_models$result)## coef 2.5 % 97.5 % vars
## (Intercept) 0.01 -0.02 0.04 (Intercept)
## scale(Pauth) 0.11 0.08 0.14 scale(Pauth)
## scale(PControl) -0.06 -0.09 -0.04 scale(PControl)
## scale(Pmon) 0.07 0.05 0.10 scale(Pmon)
## scale(PSupport) 0.18 0.15 0.21 scale(PSupport)contemporanous_models$avg_mod_summary##
## Call:
## model.avg(object = dr1)
##
## Component model call:
## lme.formula(fixed = scale(NonAttach) ~ <16 unique rhs>, data =
## ImputedDanikaData, random = ~1 | ID, correlation = corAR1(), control =
## ctrl)
##
## Component models:
## df logLik AICc delta weight
## 1234 8 -7891.94 15799.91 0.00 1
## 134 7 -7899.81 15813.64 13.74 0
## 124 7 -7904.26 15822.54 22.63 0
## 14 6 -7915.36 15842.73 42.83 0
## 234 7 -7918.23 15850.48 50.57 0
## 34 6 -7929.08 15870.17 70.27 0
## 24 6 -7939.92 15891.86 91.95 0
## 4 5 -7956.09 15922.18 122.27 0
## 123 7 -7956.27 15926.55 126.65 0
## 12 6 -7982.61 15977.24 177.33 0
## 13 6 -7993.51 15999.03 199.12 0
## 1 5 -8034.55 16079.10 279.19 0
## 23 6 -8036.67 16085.36 285.45 0
## 2 5 -8099.24 16208.49 408.59 0
## 3 5 -8106.65 16223.31 423.40 0
## (Null) 4 -8212.18 16432.36 632.45 0
##
## Term codes:
## scale(Pauth) scale(PControl) scale(Pmon) scale(PSupport)
## 1 2 3 4
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 0.01179 0.01656 0.01656 0.712 0.477
## scale(Pauth) 0.10865 0.01408 0.01408 7.715 < 0.0000000000000002
## scale(PControl) -0.06361 0.01352 0.01352 4.706 0.0000025
## scale(Pmon) 0.07324 0.01303 0.01303 5.620 < 0.0000000000000002
## scale(PSupport) 0.17850 0.01529 0.01529 11.673 < 0.0000000000000002
##
## (Intercept)
## scale(Pauth) ***
## scale(PControl) ***
## scale(Pmon) ***
## scale(PSupport) ***
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 0.01179 0.01656 0.01656 0.712 0.477
## scale(Pauth) 0.10865 0.01408 0.01408 7.715 < 0.0000000000000002
## scale(PControl) -0.06368 0.01337 0.01337 4.763 0.0000019
## scale(Pmon) 0.07324 0.01303 0.01303 5.621 < 0.0000000000000002
## scale(PSupport) 0.17850 0.01529 0.01529 11.673 < 0.0000000000000002
##
## (Intercept)
## scale(Pauth) ***
## scale(PControl) ***
## scale(Pmon) ***
## scale(PSupport) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dredge(modelRIFull)## Warning in dredge(modelRIFull): comparing models fitted by REML## Fixed term is "(Intercept)"## Global model call: lme.formula(fixed = scale(NonAttach) ~ scale(PSupport) + scale(Pmon) +
## scale(Pauth) + scale(PControl), data = ImputedDanikaData,
## random = ~1 | ID, correlation = corAR1(), control = ctrl)
## ---
## Model selection table
## (Int) scl(Pth) scl(PCn) scl(Pmn) scl(PSp) df logLik AICc delta
## 16 0.01179 0.1086 -0.06368 0.07323 0.1785 8 -7891.944 15799.9 0.00
## 14 0.01107 0.1138 0.08002 0.2025 7 -7899.814 15813.6 13.74
## 12 0.01156 0.1230 -0.07194 0.1936 7 -7904.260 15822.5 22.63
## 10 0.01093 0.1304 0.2226 6 -7915.361 15842.7 42.83
## 15 0.01191 -0.07169 0.09145 0.2234 7 -7918.233 15850.5 50.57
## 13 0.01114 0.10010 0.2528 6 -7929.081 15870.2 70.27
## 11 0.01170 -0.08365 0.2504 6 -7939.923 15891.9 91.95
## 9 0.01099 0.2883 5 -7956.085 15922.2 122.27
## 8 0.01348 0.1714 -0.11560 0.10010 7 -7956.269 15926.6 126.65
## 4 0.01361 0.1991 -0.13340 6 -7982.614 15977.2 177.33
## 6 0.01287 0.1987 0.12130 6 -7993.510 15999.0 199.12
## 2 0.01300 0.2382 5 -8034.545 16079.1 279.19
## 7 0.01459 -0.15380 0.14590 6 -8036.675 16085.4 285.45
## 3 0.01514 -0.19060 5 -8099.243 16208.5 408.59
## 5 0.01407 0.18470 5 -8106.649 16223.3 423.40
## 1 0.01497 4 -8212.176 16432.4 632.45
## weight
## 16 0.999
## 14 0.001
## 12 0.000
## 10 0.000
## 15 0.000
## 13 0.000
## 11 0.000
## 9 0.000
## 8 0.000
## 4 0.000
## 6 0.000
## 2 0.000
## 7 0.000
## 3 0.000
## 5 0.000
## 1 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | IDThis tests one random slope one at a time. probably not essential
but this shows that the estimates don’t really change.
modelRIFull <- lme(scale(NonAttach) ~ scale(PSupport)+scale(Pmon)+scale(Pauth)+ scale(PControl),
random = ~ 1 | ID,
correlation = corAR1(),
control=ctrl,
data = ImputedDanikaData)
modelRSLP1 <- update(modelRIFull, random = ~ scale(PSupport) | ID)
modelRSLP2 <- update(modelRIFull, random = ~ scale(Pmon) | ID)
modelRSLP3 <- update(modelRIFull, random = ~ scale(Pauth) | ID)
modelRSLP4 <- update(modelRIFull, random = ~ scale(PControl) | ID)
summary(modelRSLP1)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 15444.62 15513.21 -7712.309
##
## Random effects:
## Formula: ~scale(PSupport) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.5489156 (Intr)
## scale(PSupport) 0.4159799 0.193
## Residual 0.6365075
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.4619296
## Fixed effects: scale(NonAttach) ~ scale(PSupport) + scale(Pmon) + scale(Pauth) + scale(PControl)
## Value Std.Error DF t-value p-value
## (Intercept) -0.00201214 0.01642027 4690 -0.122540 0.9025
## scale(PSupport) 0.19213130 0.01835164 4690 10.469434 0.0000
## scale(Pmon) 0.05899644 0.01297361 4690 4.547420 0.0000
## scale(Pauth) 0.10257360 0.01413696 4690 7.255705 0.0000
## scale(PControl) -0.06264289 0.01325455 4690 -4.726141 0.0000
## Correlation:
## (Intr) sc(PS) scl(Pm) scl(Pt)
## scale(PSupport) 0.003
## scale(Pmon) -0.001 -0.148
## scale(Pauth) -0.006 -0.339 -0.177
## scale(PControl) -0.015 0.289 0.113 0.067
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.80495555 -0.47759089 -0.03522963 0.44203986 3.91518661
##
## Number of Observations: 7041
## Number of Groups: 2347summary(modelRSLP2)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 15482.72 15551.31 -7731.359
##
## Random effects:
## Formula: ~scale(Pmon) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.5758398 (Intr)
## scale(Pmon) 0.3989270 0.13
## Residual 0.6273648
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.4437552
## Fixed effects: scale(NonAttach) ~ scale(PSupport) + scale(Pmon) + scale(Pauth) + scale(PControl)
## Value Std.Error DF t-value p-value
## (Intercept) 0.01046827 0.01658842 4690 0.631059 0.528
## scale(PSupport) 0.17254470 0.01530421 4690 11.274327 0.000
## scale(Pmon) 0.07955029 0.01644946 4690 4.836044 0.000
## scale(Pauth) 0.10937108 0.01414480 4690 7.732245 0.000
## scale(PControl) -0.06853683 0.01336236 4690 -5.129097 0.000
## Correlation:
## (Intr) sc(PS) scl(Pm) scl(Pt)
## scale(PSupport) -0.009
## scale(Pmon) -0.019 -0.146
## scale(Pauth) -0.002 -0.377 -0.163
## scale(PControl) -0.009 0.326 0.096 0.081
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.72894923 -0.45997151 -0.03933632 0.44414380 3.77207287
##
## Number of Observations: 7041
## Number of Groups: 2347summary(modelRSLP3)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 15547.91 15616.5 -7763.955
##
## Random effects:
## Formula: ~scale(Pauth) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.5997340 (Intr)
## scale(Pauth) 0.3781102 0.051
## Residual 0.6166571
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.4082338
## Fixed effects: scale(NonAttach) ~ scale(PSupport) + scale(Pmon) + scale(Pauth) + scale(PControl)
## Value Std.Error DF t-value p-value
## (Intercept) 0.00005206 0.01656990 4690 0.003142 0.9975
## scale(PSupport) 0.17133444 0.01530823 4690 11.192307 0.0000
## scale(Pmon) 0.06122284 0.01303181 4690 4.697955 0.0000
## scale(Pauth) 0.12071941 0.01712565 4690 7.049039 0.0000
## scale(PControl) -0.05751642 0.01344147 4690 -4.279028 0.0000
## Correlation:
## (Intr) sc(PS) scl(Pm) scl(Pt)
## scale(PSupport) -0.017
## scale(Pmon) -0.001 -0.173
## scale(Pauth) -0.019 -0.328 -0.157
## scale(PControl) -0.009 0.319 0.099 0.075
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.45965163 -0.45943588 -0.03200905 0.43075788 3.99510650
##
## Number of Observations: 7041
## Number of Groups: 2347summary(modelRSLP4)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 15453.2 15521.78 -7716.598
##
## Random effects:
## Formula: ~scale(PControl) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.5550515 (Intr)
## scale(PControl) 0.4299621 -0.129
## Residual 0.6318327
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.4674289
## Fixed effects: scale(NonAttach) ~ scale(PSupport) + scale(Pmon) + scale(Pauth) + scale(PControl)
## Value Std.Error DF t-value p-value
## (Intercept) -0.00125721 0.01652264 4690 -0.076090 0.9394
## scale(PSupport) 0.18102375 0.01520552 4690 11.905137 0.0000
## scale(Pmon) 0.06488348 0.01299185 4690 4.994169 0.0000
## scale(Pauth) 0.10410889 0.01423724 4690 7.312434 0.0000
## scale(PControl) -0.07800296 0.01672137 4690 -4.664867 0.0000
## Correlation:
## (Intr) sc(PS) scl(Pm) scl(Pt)
## scale(PSupport) -0.019
## scale(Pmon) -0.004 -0.174
## scale(Pauth) -0.005 -0.377 -0.182
## scale(PControl) -0.013 0.269 0.085 0.068
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.73581553 -0.46859117 -0.02871524 0.42750612 4.29293265
##
## Number of Observations: 7041
## Number of Groups: 2347anova(modelRIFull, modelRSLP1)## Model df AIC BIC logLik Test L.Ratio p-value
## modelRIFull 1 8 15799.89 15854.76 -7891.944
## modelRSLP1 2 10 15444.62 15513.21 -7712.309 1 vs 2 359.269 <.0001anova(modelRIFull, modelRSLP2)## Model df AIC BIC logLik Test L.Ratio p-value
## modelRIFull 1 8 15799.89 15854.76 -7891.944
## modelRSLP2 2 10 15482.72 15551.31 -7731.359 1 vs 2 321.1683 <.0001anova(modelRIFull, modelRSLP3)## Model df AIC BIC logLik Test L.Ratio p-value
## modelRIFull 1 8 15799.89 15854.76 -7891.944
## modelRSLP3 2 10 15547.91 15616.50 -7763.955 1 vs 2 255.9769 <.0001anova(modelRIFull, modelRSLP4)## Model df AIC BIC logLik Test L.Ratio p-value
## modelRIFull 1 8 15799.89 15854.76 -7891.944
## modelRSLP4 2 10 15453.20 15521.78 -7716.598 1 vs 2 350.6907 <.0001Parent antecede
Lag variables
ImputedDanikaData <- ImputedDanikaData %>%
group_by(ID) %>% # data needs to be lagged within each ID
mutate(lagPmon = lag(Pmon),
lagPauth = lag(Pauth),
lagNonAttach = lag(NonAttach),
lagPControl = lag(PControl),
lagPSupport = lag(PSupport))Single model
# Then run your model:
laggedPSupportRI <- lme(scale(NonAttach) ~ scale(lagPSupport) + scale(lagNonAttach),
random = ~ 1|ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
laggedPSupportRSLP <- lme(scale(NonAttach) ~ scale(lagPSupport) + scale(lagNonAttach),
random = ~ 1+scale(lagPSupport) |ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
anova(laggedPSupportRI,laggedPSupportRSLP) #significant so we go with second model## Model df AIC BIC logLik Test L.Ratio p-value
## laggedPSupportRI 1 6 9762.714 9801.435 -4875.357
## laggedPSupportRSLP 2 8 9720.894 9772.521 -4852.447 1 vs 2 45.82027 <.0001summary(laggedPSupportRSLP)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 9720.894 9772.521 -4852.447
##
## Random effects:
## Formula: ~1 + scale(lagPSupport) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.09427226 (Intr)
## scale(lagPSupport) 0.19274682 0.677
## Residual 0.65233315
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.2694491
## Fixed effects: scale(NonAttach) ~ scale(lagPSupport) + scale(lagNonAttach)
## Value Std.Error DF t-value p-value
## (Intercept) -0.0009558 0.009025428 2346 -0.10590 0.9157
## scale(lagPSupport) 0.0397430 0.010776260 2345 3.68801 0.0002
## scale(lagNonAttach) 0.7517483 0.010066896 2345 74.67529 0.0000
## Correlation:
## (Intr) sc(PS)
## scale(lagPSupport) 0.016
## scale(lagNonAttach) 0.008 -0.340
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.411174739 -0.431381904 -0.002058009 0.424396062 7.105272218
##
## Number of Observations: 4694
## Number of Groups: 2347laggedPmonRI <- lme(scale(NonAttach) ~ scale(lagPmon) + scale(lagNonAttach),
random = ~ 1|ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
laggedPmonRSLP <- lme(scale(NonAttach) ~ scale(lagPmon) + scale(lagNonAttach),
random = ~ 1+scale(lagPmon) |ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
anova(laggedPmonRI, laggedPmonRSLP)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedPmonRI 1 6 9776.098 9814.818 -4882.049
## laggedPmonRSLP 2 8 9767.228 9818.855 -4875.614 1 vs 2 12.87009 0.0016summary(laggedPmonRSLP)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 9767.228 9818.855 -4875.614
##
## Random effects:
## Formula: ~1 + scale(lagPmon) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.03352709 (Intr)
## scale(lagPmon) 0.12303439 0.045
## Residual 0.67748640
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.2251566
## Fixed effects: scale(NonAttach) ~ scale(lagPmon) + scale(lagNonAttach)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0001004 0.009025901 2346 0.01113 0.9911
## scale(lagPmon) 0.0132725 0.010117693 2345 1.31181 0.1897
## scale(lagNonAttach) 0.7677646 0.009657924 2345 79.49582 0.0000
## Correlation:
## (Intr) scl(P)
## scale(lagPmon) -0.041
## scale(lagNonAttach) 0.000 -0.267
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.649517972 -0.443435692 -0.006868363 0.449202254 6.869131236
##
## Number of Observations: 4694
## Number of Groups: 2347laggedPauthRI <- lme(scale(NonAttach) ~ scale(lagPauth) + scale(lagNonAttach),
random = ~ 1|ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
laggedPauthRSLP <- lme(scale(NonAttach) ~ scale(lagPauth) + scale(lagNonAttach),
random = ~ 1+scale(lagPauth) |ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
anova(laggedPauthRI, laggedPauthRSLP)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedPauthRI 1 6 9771.280 9810.000 -4879.640
## laggedPauthRSLP 2 8 9762.976 9814.604 -4873.488 1 vs 2 12.30366 0.0021summary(laggedPauthRSLP)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 9762.976 9814.604 -4873.488
##
## Random effects:
## Formula: ~1 + scale(lagPauth) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.0406079 (Intr)
## scale(lagPauth) 0.1334206 0.329
## Residual 0.6748081
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.227026
## Fixed effects: scale(NonAttach) ~ scale(lagPauth) + scale(lagNonAttach)
## Value Std.Error DF t-value p-value
## (Intercept) -0.0021221 0.009028588 2346 -0.23504 0.8142
## scale(lagPauth) 0.0275668 0.010475771 2345 2.63148 0.0086
## scale(lagNonAttach) 0.7611387 0.009902916 2345 76.86005 0.0000
## Correlation:
## (Intr) scl(P)
## scale(lagPauth) -0.020
## scale(lagNonAttach) 0.009 -0.326
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.626482487 -0.451393190 -0.007700151 0.449238297 6.960189188
##
## Number of Observations: 4694
## Number of Groups: 2347laggedPControlRI <- lme(scale(NonAttach) ~ scale(lagPControl) + scale(lagNonAttach),
random = ~ 1|ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
laggedPControlRSLP <- lme(scale(NonAttach) ~ scale(lagPControl) + scale(lagNonAttach),
random = ~ 1+scale(lagPControl) |ID,
correlation = corAR1(),
na.action=na.exclude,
control=ctrl,
data = ImputedDanikaData)
anova(laggedPControlRI, laggedPControlRSLP)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedPControlRI 1 6 9776.858 9815.578 -4882.429
## laggedPControlRSLP 2 8 9737.171 9788.799 -4860.586 1 vs 2 43.6865 <.0001summary(laggedPControlRSLP)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 9737.171 9788.799 -4860.586
##
## Random effects:
## Formula: ~1 + scale(lagPControl) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.1453872 (Intr)
## scale(lagPControl) 0.1848083 -0.048
## Residual 0.6468418
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.3108885
## Fixed effects: scale(NonAttach) ~ scale(lagPControl) + scale(lagNonAttach)
## Value Std.Error DF t-value p-value
## (Intercept) -0.0027568 0.009012519 2346 -0.30589 0.7597
## scale(lagPControl) -0.0134994 0.010641637 2345 -1.26855 0.2047
## scale(lagNonAttach) 0.7672460 0.009815915 2345 78.16347 0.0000
## Correlation:
## (Intr) sc(PC)
## scale(lagPControl) 0.052
## scale(lagNonAttach) 0.007 0.297
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.864824000 -0.416932006 0.003372449 0.415270322 7.206096098
##
## Number of Observations: 4694
## Number of Groups: 2347Graph it
# Extract fixed effects
fixed_effects <- fixef(laggedPSupportRSLP)
# Extract random effects
random_effects <- ranef(laggedPSupportRSLP)
# 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(lagPSupport)`
# Add the fixed effect to each random slope to get the slope for each person
random_effects_df$person_slope <- fixed_effects[2] + random_effects_df$random_slope
# Calculate the mean of person_slope
mean_person_slope2 <- mean(random_effects_df$person_slope, na.rm = TRUE)
# Create the histogram
Laghistogram <- 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_slope2), color="red", linetype="dashed", size=1) +
annotate("text", x = mean_person_slope2, y= -Inf, label = paste("Mean =", round(mean_person_slope2, 2)), vjust = -1.5) +
theme_minimal() +
labs(x="Person Slope", y="Count", title="Within person association: PSupport with Nonattachment")
# Print the histogram
print(Laghistogram)
All variables at once
need to remove nas at start of lagged time series to run information theoretic model,ing
ImputedDanikaData2 <- ImputedDanikaData %>%
group_by(ID) %>% # Data needs to be lagged within each ID
mutate(lagPmon = lag(Pmon),
lagPauth = lag(Pauth),
lagNonAttach = lag(NonAttach),
lagPControl = lag(PControl),
lagPSupport = lag(PSupport)) %>%
drop_na(lagPmon, lagPauth, lagNonAttach, lagPControl, lagPSupport) # Drop rows with NA in these columns
laggedTotalRSLP <- lme(
scale(NonAttach) ~ scale(lagNonAttach) + scale(lagPSupport) + scale(lagPControl) + scale(lagPauth)+
lag(Pmon),
random = ~ 1 | ID,
correlation = corAR1(),
control = ctrl,
data = ImputedDanikaData2
)identify the most important lagged variable using information theoretic modelling
lagged_models<- information_theoretic_tables(laggedTotalRSLP)## Warning in dredge(model): comparing models fitted by REML## Fixed term is "(Intercept)"print(lagged_models$var_importance)## scale(lagNonAttach) scale(lagPSupport) scale(lagPauth)
## Sum of weights: 1 0.93 0.02
## N containing models: 16 16 16
## scale(lagPControl) lag(Pmon)
## Sum of weights: 0.02 <0.01
## N containing models: 16 16print(lagged_models$model_selection)## Global model call: lme.formula(fixed = scale(NonAttach) ~ scale(lagNonAttach) +
## scale(lagPSupport) + scale(lagPControl) + scale(lagPauth) +
## lag(Pmon), data = ImputedDanikaData2, random = ~1 | ID, correlation = corAR1(),
## control = ctrl)
## ---
## Model selection table
## (Int) lag(Pmn) scl(lNA) scl(lgP) scl(lPC) scl(lPS) df
## 19 0.00000000000000023460 0.75910 0.03840 6
## 3 0.00000000000000025880 0.77410 5
## 27 0.00000000000000024490 0.76140 0.014340 0.04605 7
## 7 0.00000000000000030360 0.76460 0.025120 6
## 23 0.00000000000000024710 0.75840 0.003995 0.03604 7
## 11 0.00000000000000029090 0.77040 -0.009744 6
## 31 0.00000000000000026100 0.76070 0.004868 0.014590 0.04330 8
## 15 0.00000000000000024070 0.76430 0.024750 -0.001089 7
## 4 0.04920999999999999680 -0.019590 0.72390 6
## 29 0.00000000000000004776 0.092770 -0.076380 0.09612 7
## 21 0.00000000000000009376 0.099400 0.13020 6
## 20 0.05378999999999999754 -0.021580 0.06060 0.18820 7
## 13 0.00000000000000010030 0.132300 -0.109000 6
## 25 0.00000000000000005541 -0.083960 0.14320 6
## 17 -0.00000000000000002226 0.18440 5
## 8 0.05802000000000000213 -0.023340 0.05596 0.167300 7
## 26 0.07631000000000000283 -0.030870 -0.003955 0.18400 7
## 22 0.07671000000000000041 -0.031050 0.003414 0.18410 7
## 5 -0.00000000000000002418 0.166100 5
## 12 0.05414000000000000062 -0.021750 0.05570 -0.151600 7
## 32 0.01626000000000000015 -0.005848 -0.01301 -0.001881 -0.004351 0.18370 9
## 9 0.00000000000000017380 -0.150100 5
## 14 0.07198000000000000231 -0.029130 0.006293 -0.149500 7
## 1 0.00000000000000002097 4
## 2 -0.03257000000000000173 0.013770 5
## 18 0.01532000000000000028 -0.005695 0.30410 6
## 24 0.03540000000000000091 -0.013910 -0.01195 0.027310 0.29660 8
## 30 0.03504000000000000170 -0.013730 -0.014300 -0.024490 0.29800 8
## 28 0.03481000000000000066 -0.013640 -0.01123 -0.023440 0.29810 8
## 6 0.06393999999999999684 -0.025830 0.272200 6
## 16 0.03885999999999999871 -0.015310 -0.01545 0.063040 -0.237200 8
## 10 0.01992999999999999980 -0.007829 -0.249100 6
## logLik AICc delta weight
## 19 -4875.357 9762.7 0.00 0.903
## 3 -4879.211 9768.4 5.70 0.052
## 27 -4878.129 9770.3 7.55 0.021
## 7 -4879.640 9771.3 8.57 0.012
## 23 -4878.797 9771.6 8.89 0.011
## 11 -4882.429 9776.9 14.14 0.001
## 31 -4881.542 9779.1 16.38 0.000
## 15 -4883.281 9780.6 17.85 0.000
## 4 -4924.949 9861.9 99.18 0.000
## 29 -5739.576 11493.2 1730.44 0.000
## 21 -5746.065 11504.1 1741.42 0.000
## 20 -5746.706 11507.4 1744.70 0.000
## 13 -5748.164 11508.3 1745.61 0.000
## 25 -5749.775 11511.6 1748.84 0.000
## 17 -5758.274 11526.6 1763.83 0.000
## 8 -5758.543 11531.1 1768.38 0.000
## 26 -5761.837 11537.7 1774.97 0.000
## 22 -5761.852 11537.7 1775.00 0.000
## 5 -5767.758 11545.5 1782.80 0.000
## 12 -5769.462 11552.9 1790.22 0.000
## 32 -5769.299 11556.6 1793.90 0.000
## 9 -5778.574 11567.2 1804.43 0.000
## 14 -5782.326 11578.7 1815.94 0.000
## 1 -5821.443 11650.9 1888.16 0.000
## 2 -5824.506 11659.0 1896.29 0.000
## 18 -6439.680 12891.4 3128.65 0.000
## 24 -6444.002 12904.0 3141.30 0.000
## 30 -6444.207 12904.4 3141.71 0.000
## 28 -6444.462 12905.0 3142.22 0.000
## 6 -6486.555 12985.1 3222.40 0.000
## 16 -6513.610 13043.3 3280.52 0.000
## 10 -6517.065 13046.1 3283.42 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## 1 | IDprint(lagged_models$result)## coef 2.5 % 97.5 % vars
## (Intercept) 0.00 -0.02 0.02 (Intercept)
## scale(lagNonAttach) 0.76 0.74 0.78 scale(lagNonAttach)
## scale(lagPSupport) 0.04 0.01 0.06 scale(lagPSupport)
## scale(lagPControl) 0.00 0.00 0.01 scale(lagPControl)
## scale(lagPauth) 0.00 -0.01 0.01 scale(lagPauth)
## lag(Pmon) 0.00 0.00 0.00 lag(Pmon)lagged_models$avg_mod_summary##
## Call:
## model.avg(object = dr1)
##
## Component model call:
## lme.formula(fixed = scale(NonAttach) ~ <32 unique rhs>, data =
## ImputedDanikaData2, random = ~1 | ID, correlation = corAR1(), control =
## ctrl)
##
## Component models:
## df logLik AICc delta weight
## 25 6 -4875.36 9762.73 0.00 0.90
## 2 5 -4879.21 9768.44 5.70 0.05
## 245 7 -4878.13 9770.28 7.55 0.02
## 23 6 -4879.64 9771.30 8.57 0.01
## 235 7 -4878.80 9771.62 8.89 0.01
## 24 6 -4882.43 9776.88 14.14 0.00
## 2345 8 -4881.54 9779.11 16.38 0.00
## 234 7 -4883.28 9780.59 17.85 0.00
## 12 6 -4924.95 9861.92 99.18 0.00
## 345 7 -5739.58 11493.18 1730.44 0.00
## 35 6 -5746.06 11504.15 1741.42 0.00
## 125 7 -5746.71 11507.44 1744.70 0.00
## 34 6 -5748.16 11508.35 1745.61 0.00
## 45 6 -5749.77 11511.57 1748.84 0.00
## 5 5 -5758.27 11526.56 1763.83 0.00
## 123 7 -5758.54 11531.11 1768.38 0.00
## 145 7 -5761.84 11537.70 1774.97 0.00
## 135 7 -5761.85 11537.73 1775.00 0.00
## 3 5 -5767.76 11545.53 1782.80 0.00
## 124 7 -5769.46 11552.95 1790.22 0.00
## 12345 9 -5769.30 11556.64 1793.90 0.00
## 4 5 -5778.57 11567.16 1804.43 0.00
## 134 7 -5782.33 11578.68 1815.94 0.00
## (Null) 4 -5821.44 11650.89 1888.16 0.00
## 1 5 -5824.51 11659.02 1896.29 0.00
## 15 6 -6439.68 12891.38 3128.65 0.00
## 1235 8 -6444.00 12904.03 3141.30 0.00
## 1345 8 -6444.21 12904.44 3141.71 0.00
## 1245 8 -6444.46 12904.95 3142.22 0.00
## 13 6 -6486.56 12985.13 3222.40 0.00
## 1234 8 -6513.61 13043.25 3280.52 0.00
## 14 6 -6517.06 13046.15 3283.42 0.00
##
## Term codes:
## lag(Pmon) scale(lagNonAttach) scale(lagPauth) scale(lagPControl)
## 1 2 3 4
## scale(lagPSupport)
## 5
##
## Model-averaged coefficients:
## (full average)
## Estimate
## (Intercept) 0.000000000000000237159081741
## scale(lagNonAttach) 0.759976207552177651649571999
## scale(lagPSupport) 0.036021387255877078403809577
## scale(lagPControl) 0.000293032570565425021062722
## scale(lagPauth) 0.000359702887321991320555836
## lag(Pmon) -0.000000000000000000000005131
## Std. Error
## (Intercept) 0.008998500851997786112446143
## scale(lagNonAttach) 0.010448334789386773055475111
## scale(lagPSupport) 0.013583173454400349414172844
## scale(lagPControl) 0.002669405720011999803148450
## scale(lagPauth) 0.003284390884447106577492814
## lag(Pmon) 0.000000000000438443036467207
## Adjusted SE z value Pr(>|z|)
## (Intercept) 0.009003145776820606696611016 0.000 1.00000
## scale(lagNonAttach) 0.010453159629118071935471868 72.703 < 0.0000000000000002
## scale(lagPSupport) 0.013586683484904186688946304 2.651 0.00802
## scale(lagPControl) 0.002669953028366566991080866 0.110 0.91261
## scale(lagPauth) 0.003284831563426439393627021 0.110 0.91280
## lag(Pmon) 0.000000000000438551102506773 0.000 1.00000
##
## (Intercept)
## scale(lagNonAttach) ***
## scale(lagPSupport) **
## scale(lagPControl)
## scale(lagPauth)
## lag(Pmon)
##
## (conditional average)
## Estimate Std. Error
## (Intercept) 0.0000000000000002372 0.0089985008519977861
## scale(lagNonAttach) 0.7599762075521776516 0.0104483347893867731
## scale(lagPSupport) 0.0385458525976607250 0.0100062639180349446
## scale(lagPControl) 0.0134111660849120587 0.0122552978053976011
## scale(lagPauth) 0.0153388521971933086 0.0151735912814050279
## lag(Pmon) -0.0195933242685443695 0.0187126592021831795
## Adjusted SE z value Pr(>|z|)
## (Intercept) 0.0090031457768206067 0.000 1.000000
## scale(lagNonAttach) 0.0104531596291180719 72.703 < 0.0000000000000002 ***
## scale(lagPSupport) 0.0100113619540154099 3.850 0.000118 ***
## scale(lagPControl) 0.0122607531389596347 1.094 0.274030
## scale(lagPauth) 0.0151776586043956257 1.011 0.312198
## lag(Pmon) 0.0187223267120488895 1.047 0.295320
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We already established support significantly varies when analyzied by itself. here we see no eviendence of anything else being random, but these are variables together
model1 <- update(laggedTotalRSLP, random = ~ scale(lagNonAttach) | ID)
model2 <- update(laggedTotalRSLP, random = ~ scale(lagPSupport) | ID)
model3 <- update(laggedTotalRSLP, random = ~ scale(lagPControl) | ID)
model4 <- update(laggedTotalRSLP, random = ~ scale(lagPauth) | ID)
anova(laggedTotalRSLP, model1)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedTotalRSLP 1 9 11685.13 11743.21 -5833.567
## model1 2 11 11686.68 11757.66 -5832.339 1 vs 2 2.454814 0.2931anova(laggedTotalRSLP, model2)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedTotalRSLP 1 9 11685.13 11743.21 -5833.567
## model2 2 11 11689.25 11760.23 -5833.624 1 vs 2 0.1142883 0.9445anova(laggedTotalRSLP, model3)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedTotalRSLP 1 9 11685.13 11743.21 -5833.567
## model3 2 11 11689.27 11760.25 -5833.636 1 vs 2 0.1387022 0.933anova(laggedTotalRSLP, model4)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedTotalRSLP 1 9 11685.13 11743.21 -5833.567
## model4 2 11 11689.18 11760.16 -5833.590 1 vs 2 0.04748387 0.9765final model with only support being important
finalLagRSLP <- lme(
scale(NonAttach) ~ scale(lagNonAttach) + scale(lagPSupport),
random = ~ 1+scale(lagPSupport) | ID,
correlation = corAR1(),
control = ctrl,
data = ImputedDanikaData2
)
summary(finalLagRSLP)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData2
## AIC BIC logLik
## 9720.894 9772.521 -4852.447
##
## Random effects:
## Formula: ~1 + scale(lagPSupport) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.09427226 (Intr)
## scale(lagPSupport) 0.19274682 0.677
## Residual 0.65233315
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.2694491
## Fixed effects: scale(NonAttach) ~ scale(lagNonAttach) + scale(lagPSupport)
## Value Std.Error DF t-value p-value
## (Intercept) -0.0009558 0.009025428 2346 -0.10590 0.9157
## scale(lagNonAttach) 0.7517483 0.010066896 2345 74.67529 0.0000
## scale(lagPSupport) 0.0397430 0.010776260 2345 3.68801 0.0002
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) 0.008
## scale(lagPSupport) 0.016 -0.340
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.411174739 -0.431381904 -0.002058009 0.424396062 7.105272218
##
## Number of Observations: 4694
## Number of Groups: 2347parent conseq
one variable at time
looks like nonattached kids like there parents style of parenting more in the next year
note each dv is different in this analysis so we are not pitting each variable against the other and not using information theoretic modeling to identify most importat
laggedPmonModel <- lme(scale(Pmon) ~ scale(lagNonAttach) + scale(lagPmon),
random = ~ 1 | ID,
correlation = corAR1(),
na.action = na.exclude,
control = ctrl,
data = ImputedDanikaData)
laggedPauthModel <- lme(scale(Pauth) ~ scale(lagNonAttach) + scale(lagPauth),
random = ~ 1 | ID,
correlation = corAR1(),
na.action = na.exclude,
control = ctrl,
data = ImputedDanikaData)
laggedPControlModel <- lme(scale(PControl) ~ scale(lagNonAttach) + scale(lagPControl),
random = ~ 1 | ID,
correlation = corAR1(),
na.action = na.exclude,
control = ctrl,
data = ImputedDanikaData)
laggedPSupportModel <- lme(scale(PSupport) ~ scale(lagNonAttach) + scale(lagPSupport),
random = ~ 1 | ID,
correlation = corAR1(),
na.action = na.exclude,
control = ctrl,
data = ImputedDanikaData)
summary(laggedPmonModel)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 9029.693 9068.413 -4508.846
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.0777084 0.6285826
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.1398631
## Fixed effects: scale(Pmon) ~ scale(lagNonAttach) + scale(lagPmon)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0000000 0.008658792 2346 0.00000 1.0000
## scale(lagNonAttach) 0.0344614 0.009166749 2345 3.75940 0.0002
## scale(lagPmon) 0.7896881 0.009148687 2345 86.31709 0.0000
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) 0.000
## scale(lagPmon) 0.000 -0.284
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.6834341 -0.3668264 0.1145957 0.3652129 5.4730292
##
## Number of Observations: 4694
## Number of Groups: 2347summary(laggedPauthModel)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 8980.326 9019.047 -4484.163
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.1963383 0.6007247
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.2893636
## Fixed effects: scale(Pauth) ~ scale(lagNonAttach) + scale(lagPauth)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0000000 0.008429574 2346 0.00000 1
## scale(lagNonAttach) 0.0429264 0.009140190 2345 4.69645 0
## scale(lagPauth) 0.7949602 0.009115923 2345 87.20567 0
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) 0.00
## scale(lagPauth) 0.00 -0.34
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.01932036 -0.38169837 0.03460203 0.37727319 5.31125572
##
## Number of Observations: 4694
## Number of Groups: 2347summary(laggedPControlModel)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 9076.257 9114.977 -4532.128
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.05088582 0.6371803
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.1852208
## Fixed effects: scale(PControl) ~ scale(lagNonAttach) + scale(lagPControl)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0000000 0.008460264 2346 0.00000 1
## scale(lagNonAttach) -0.0401622 0.009140306 2345 -4.39396 0
## scale(lagPControl) 0.7914212 0.009107913 2345 86.89380 0
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) 0.000
## scale(lagPControl) 0.000 0.327
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.75589873 -0.41564841 -0.05596209 0.39269471 5.81544683
##
## Number of Observations: 4694
## Number of Groups: 2347summary(laggedPSupportModel)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 8662.105 8700.825 -4325.052
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02891586 0.614166
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.2304892
## Fixed effects: scale(PSupport) ~ scale(lagNonAttach) + scale(lagPSupport)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0000000 0.007886223 2346 0.00000 1.0000
## scale(lagNonAttach) 0.0189736 0.008714016 2345 2.17737 0.0296
## scale(lagPSupport) 0.8272891 0.008666196 2345 95.46163 0.0000
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) 0.000
## scale(lagPSupport) 0.000 -0.371
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.27229688 -0.41116549 0.08809055 0.40355502 5.72570350
##
## Number of Observations: 4694
## Number of Groups: 2347Testing for random slopes
all random except for monitoring
# Random Slope Model for Pmon
laggedPmonModelRS <- lme(scale(Pmon) ~ scale(lagNonAttach) + scale(lagPmon),
random = ~ 1+scale(lagNonAttach) | ID,
correlation = corAR1(),
na.action = na.exclude,
control = ctrl,
data = ImputedDanikaData)
# Random Slope Model for Pauth
laggedPauthModelRS <- lme(scale(Pauth) ~ scale(lagNonAttach) + scale(lagPauth),
random = ~ 1+scale(lagNonAttach) | ID,
correlation = corAR1(),
na.action = na.exclude,
control = ctrl,
data = ImputedDanikaData)
# Random Slope Model for PControl
laggedPControlModelRS <- lme(scale(PControl) ~ scale(lagNonAttach) + scale(lagPControl),
random = ~ 1+scale(lagNonAttach) | ID,
correlation = corAR1(),
na.action = na.exclude,
control = ctrl,
data = ImputedDanikaData)
# Random Slope Model for PSupport
laggedPSupportModelRS <- lme(scale(PSupport) ~ scale(lagNonAttach) + scale(lagPSupport),
random = ~ 1+scale(lagNonAttach) | ID,
correlation = corAR1(),
na.action = na.exclude,
control = ctrl,
data = ImputedDanikaData)
# ANOVA test for Pmon
anova(laggedPmonModel, laggedPmonModelRS)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedPmonModel 1 6 9029.693 9068.413 -4508.846
## laggedPmonModelRS 2 8 8982.995 9034.622 -4483.498 1 vs 2 50.69752 <.0001# ANOVA test for Pauth
anova(laggedPauthModel, laggedPauthModelRS)## Model df AIC BIC logLik Test L.Ratio p-value
## laggedPauthModel 1 6 8980.326 9019.047 -4484.163
## laggedPauthModelRS 2 8 8964.891 9016.518 -4474.445 1 vs 2 19.43536 0.0001# ANOVA test for PControl
anova(laggedPControlModel, laggedPControlModelRS)## Model df AIC BIC logLik Test L.Ratio
## laggedPControlModel 1 6 9076.257 9114.977 -4532.128
## laggedPControlModelRS 2 8 9075.886 9127.513 -4529.943 1 vs 2 4.371148
## p-value
## laggedPControlModel
## laggedPControlModelRS 0.1124# ANOVA test for PSupport
anova(laggedPSupportModel, laggedPSupportModelRS)## Model df AIC BIC logLik Test L.Ratio
## laggedPSupportModel 1 6 8662.105 8700.825 -4325.052
## laggedPSupportModelRS 2 8 8607.848 8659.476 -4295.924 1 vs 2 58.25624
## p-value
## laggedPSupportModel
## laggedPSupportModelRS <.0001consequential models with random slopes except for parental monitoring , whcih was not significant
summary(laggedPmonModel)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 9029.693 9068.413 -4508.846
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.0777084 0.6285826
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.1398631
## Fixed effects: scale(Pmon) ~ scale(lagNonAttach) + scale(lagPmon)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0000000 0.008658792 2346 0.00000 1.0000
## scale(lagNonAttach) 0.0344614 0.009166749 2345 3.75940 0.0002
## scale(lagPmon) 0.7896881 0.009148687 2345 86.31709 0.0000
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) 0.000
## scale(lagPmon) 0.000 -0.284
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.6834341 -0.3668264 0.1145957 0.3652129 5.4730292
##
## Number of Observations: 4694
## Number of Groups: 2347summary(laggedPauthModelRS)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 8964.891 9016.518 -4474.445
##
## Random effects:
## Formula: ~1 + scale(lagNonAttach) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.1344928 (Intr)
## scale(lagNonAttach) 0.1215738 -0.5
## Residual 0.6049685
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.2452998
## Fixed effects: scale(Pauth) ~ scale(lagNonAttach) + scale(lagPauth)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0006413 0.008442787 2346 0.07595 0.9395
## scale(lagNonAttach) 0.0431825 0.009691282 2345 4.45581 0.0000
## scale(lagPauth) 0.7926007 0.009177487 2345 86.36358 0.0000
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) -0.072
## scale(lagPauth) -0.002 -0.322
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.17187847 -0.39118329 0.03463666 0.39175288 4.85279239
##
## Number of Observations: 4694
## Number of Groups: 2347summary(laggedPControlModelRS)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 9075.886 9127.513 -4529.943
##
## Random effects:
## Formula: ~1 + scale(lagNonAttach) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.19190012 (Intr)
## scale(lagNonAttach) 0.05636297 -0.44
## Residual 0.60709379
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.3032404
## Fixed effects: scale(PControl) ~ scale(lagNonAttach) + scale(lagPControl)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0000734 0.008457370 2346 0.00867 0.9931
## scale(lagNonAttach) -0.0417583 0.009271961 2345 -4.50372 0.0000
## scale(lagPControl) 0.7912657 0.009142283 2345 86.55012 0.0000
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) -0.048
## scale(lagPControl) 0.001 0.325
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -5.15104407 -0.37651205 -0.04299083 0.35375338 5.78547918
##
## Number of Observations: 4694
## Number of Groups: 2347summary(laggedPSupportModelRS)## Linear mixed-effects model fit by REML
## Data: ImputedDanikaData
## AIC BIC logLik
## 8607.848 8659.476 -4295.924
##
## Random effects:
## Formula: ~1 + scale(lagNonAttach) | ID
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.1579808 (Intr)
## scale(lagNonAttach) 0.1163936 -0.836
## Residual 0.5818131
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## -0.3487982
## Fixed effects: scale(PSupport) ~ scale(lagNonAttach) + scale(lagPSupport)
## Value Std.Error DF t-value p-value
## (Intercept) 0.0012862 0.007866212 2346 0.16351 0.8701
## scale(lagNonAttach) 0.0159769 0.009071651 2345 1.76119 0.0783
## scale(lagPSupport) 0.8309053 0.008706855 2345 95.43117 0.0000
## Correlation:
## (Intr) sc(NA)
## scale(lagNonAttach) -0.155
## scale(lagPSupport) -0.004 -0.354
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.45505113 -0.38175310 0.07484786 0.38400862 4.98709537
##
## Number of Observations: 4694
## Number of Groups: 2347nonattach does not significant link to later parenting, as fixed effect, but there is significant heterogenity. lets look at that
# Convert the random effects to a data frame
random_effects_df <- as.data.frame(ranef(laggedPSupportModelRS))
# Isolate the random effects for scale(lagNonAttach)
random_effects_df_lagNonAttach <- random_effects_df$`scale(lagNonAttach)`
# Create the histogram
histogram <- ggplot(data.frame(lagNonAttach = random_effects_df_lagNonAttach), aes(x = lagNonAttach)) +
geom_histogram(color="black", fill="lightblue", bins=30) +
theme_minimal() +
labs(x="Random Effect of lagNonAttach", y="Count", title="Histogram of Random Effects for lagNonAttach")
# Print the histogram
print(histogram)