# import data
library(readr)
d <- read_csv("/Users/melissalagunas/Desktop/Lab/Contextualizing\ Stigma\ and\ Trauma/LatineData_Subscales.csv")
## New names:
## Rows: 259 Columns: 16
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," dbl
## (16): ...1, SSRPH, SSOSH, MHSAS, PR_PF, PR_A, PR_D, PR_SR, PR_PI, ANG, C...
## ℹ 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.
## • `` -> `...1`
View(d)
# import data
library(readr)
dat <- read_csv("/Users/melissalagunas/Desktop/Lab/Contextualizing\ Stigma\ and\ Trauma/LatineData2.csv")
## New names:
## Rows: 259 Columns: 151
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (4): AHELP, GEN_8_TEXT, ETH_7_TEXT, BORN_1 dbl (146): ...1, Progress,
## Duration..in.seconds., Consent, SSRPH_1, SSRPH_2,... lgl (1): EDU_11_TEXT
## ℹ 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.
## • `` -> `...1`
View(dat)
# Extract variables of interest
dat <- dat[, 88:100]
# Recode
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ purrr 1.0.2
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Recode the scale with dichotomous items
dat$scale <- rowSums(dat[, c("SLESQ_1", "SLESQ_2", "SLESQ_3", "SLESQ_4", "SLESQ_5", "SLESQ_6", "SLESQ_7", "SLESQ_8", "SLESQ_9", "SLESQ_10", "SLESQ_11","SLESQ_12", "SLESQ_13")])
#data set with only familismo (PHFS), stressful life events (SLESQ), and central sensitization (CSA)
library(dplyr)
d <- data.frame(d$PHFS, dat$scale, d$CSA, d$PCL5)
d <- rename(d, familism = d.PHFS, stressful_life_events = dat.scale, central_sensitization = d.CSA, trauma_symptoms = d.PCL5 )
# remove NA variables
d <- na.omit(d)
# loading packages
library(lavaan)
## This is lavaan 0.6-16
## lavaan is FREE software! Please report any bugs.
library(semPlot)
library(tidyverse)
library(psych)
##
## Attaching package: 'psych'
## The following object is masked from 'package:lavaan':
##
## cor2cov
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(jtools)
library(interactions)
library(ggplot2)
library(tidySEM)
## Loading required package: OpenMx
## OpenMx may run faster if it is compiled to take advantage of multiple cores.
##
## Attaching package: 'OpenMx'
## The following object is masked from 'package:psych':
##
## tr
## Registered S3 method overwritten by 'tidySEM':
## method from
## predict.MxModel OpenMx
##
## Attaching package: 'tidySEM'
## The following object is masked from 'package:jtools':
##
## get_data
# Correlation analysis
table <- apaTables::apa.cor.table(d, table.number = 1, show.sig.stars = TRUE,
landscape = TRUE, filename = "table.doc")
## Registered S3 methods overwritten by 'broom':
## method from
## tidy.glht jtools
## tidy.summary.glht jtools
print(table)
##
##
## Table 1
##
## Means, standard deviations, and correlations with confidence intervals
##
##
## Variable M SD 1 2 3
## 1. familism 3.67 0.96
##
## 2. stressful_life_events 16.72 2.85 -.21**
## [-.32, -.09]
##
## 3. central_sensitization 1.63 0.72 -.34** .33**
## [-.45, -.23] [.21, .43]
##
## 4. trauma_symptoms 1.51 0.97 -.37** .30** .68**
## [-.47, -.25] [.18, .40] [.61, .75]
##
##
## Note. M and SD are used to represent mean and standard deviation, respectively.
## Values in square brackets indicate the 95% confidence interval.
## The confidence interval is a plausible range of population correlations
## that could have caused the sample correlation (Cumming, 2014).
## * indicates p < .05. ** indicates p < .01.
##
psych::pairs.panels(d)
# Checking the table of means, standard deviations, and correlations
table <- apaTables::apa.cor.table(d, table.number = 1, show.sig.stars = TRUE,
landscape = TRUE, filename = "Lewis_Corr.doc")
print(table)
##
##
## Table 1
##
## Means, standard deviations, and correlations with confidence intervals
##
##
## Variable M SD 1 2 3
## 1. familism 3.67 0.96
##
## 2. stressful_life_events 16.72 2.85 -.21**
## [-.32, -.09]
##
## 3. central_sensitization 1.63 0.72 -.34** .33**
## [-.45, -.23] [.21, .43]
##
## 4. trauma_symptoms 1.51 0.97 -.37** .30** .68**
## [-.47, -.25] [.18, .40] [.61, .75]
##
##
## Note. M and SD are used to represent mean and standard deviation, respectively.
## Values in square brackets indicate the 95% confidence interval.
## The confidence interval is a plausible range of population correlations
## that could have caused the sample correlation (Cumming, 2014).
## * indicates p < .05. ** indicates p < .01.
##
# quick peak at data
psych::describe(d)
## vars n mean sd median trimmed mad min max range
## familism 1 251 3.67 0.96 3.80 3.76 0.89 1 5.0 4.0
## stressful_life_events 2 251 16.72 2.85 16.00 16.48 2.97 13 26.0 13.0
## central_sensitization 3 251 1.63 0.72 1.64 1.63 0.77 0 3.8 3.8
## trauma_symptoms 4 251 1.51 0.97 1.50 1.48 1.11 0 4.0 4.0
## skew kurtosis se
## familism -0.74 0.18 0.06
## stressful_life_events 0.69 -0.03 0.18
## central_sensitization 0.03 -0.31 0.05
## trauma_symptoms 0.26 -0.83 0.06
# Y = central sensitization; X = stressful life events; W = familism
Mod_a_path <- lm(central_sensitization ~ stressful_life_events* familism, data = d)
jtools::summ(Mod_a_path, digits = 3)
## MODEL INFO:
## Observations: 251
## Dependent Variable: central_sensitization
## Type: OLS linear regression
##
## MODEL FIT:
## F(3,247) = 18.539, p = 0.000
## R² = 0.184
## Adj. R² = 0.174
##
## Standard errors: OLS
## ----------------------------------------------------------------------
## Est. S.E. t val. p
## ------------------------------------ -------- ------- -------- -------
## (Intercept) 1.120 1.102 1.016 0.310
## stressful_life_events 0.077 0.062 1.235 0.218
## familism -0.169 0.284 -0.595 0.553
## stressful_life_events:familism -0.003 0.016 -0.161 0.872
## ----------------------------------------------------------------------
#graph interaction
interactions::interact_plot(Mod_a_path, pred = stressful_life_events, modx = familism) +
ylim(0, 10)
#probing the interaction with simple slopes
interactions::sim_slopes(Mod_a_path, pred = stressful_life_events, modx = familism)
## JOHNSON-NEYMAN INTERVAL
##
## When familism is INSIDE the interval [1.64, 5.45], the slope of
## stressful_life_events is p < .05.
##
## Note: The range of observed values of familism is [1.00, 5.00]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of stressful_life_events when familism = 2.709182 (- 1 SD):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.07 0.02 3.15 0.00
##
## Slope of stressful_life_events when familism = 3.672510 (Mean):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.07 0.01 4.52 0.00
##
## Slope of stressful_life_events when familism = 4.635838 (+ 1 SD):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.06 0.02 3.08 0.00
# Y = trauma symtpoms; X = stressful life events; W = familism
Mod_b_path <- lm(trauma_symptoms ~ stressful_life_events * familism, data = d)
jtools::summ(Mod_b_path, digits = 3)
## MODEL INFO:
## Observations: 251
## Dependent Variable: trauma_symptoms
## Type: OLS linear regression
##
## MODEL FIT:
## F(3,247) = 19.436, p = 0.000
## R² = 0.191
## Adj. R² = 0.181
##
## Standard errors: OLS
## ----------------------------------------------------------------------
## Est. S.E. t val. p
## ------------------------------------ -------- ------- -------- -------
## (Intercept) 3.290 1.479 2.225 0.027
## stressful_life_events -0.031 0.084 -0.375 0.708
## familism -0.830 0.381 -2.179 0.030
## stressful_life_events:familism 0.029 0.022 1.348 0.179
## ----------------------------------------------------------------------
#graph interaction
interactions::interact_plot(Mod_b_path, pred = stressful_life_events, modx = familism) +
ylim(1, 10)
## Warning: Removed 14 rows containing missing values (`geom_path()`).
#probing the interaction with simple slopes
interactions::sim_slopes(Mod_b_path, pred = stressful_life_events, modx = familism)
## JOHNSON-NEYMAN INTERVAL
##
## When familism is INSIDE the interval [2.88, 9.24], the slope of
## stressful_life_events is p < .05.
##
## Note: The range of observed values of familism is [1.00, 5.00]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of stressful_life_events when familism = 2.709182 (- 1 SD):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.05 0.03 1.63 0.11
##
## Slope of stressful_life_events when familism = 3.672510 (Mean):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.08 0.02 3.84 0.00
##
## Slope of stressful_life_events when familism = 4.635838 (+ 1 SD):
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.11 0.03 3.72 0.00
# Y = central sensitization; X = stressful life events; M = trauma symptoms
LMedModel <- "
central_sensitization ~ b*trauma_symptoms + c_p*stressful_life_events
trauma_symptoms ~a*stressful_life_events
#intercepts
trauma_symptoms ~ trauma_symptoms.mean*1
central_sensitization ~ central_sensitization.mean*1
indirect := a*b
direct := c_p
total_c := c_p + (a*b)
"
set.seed(230925) #required for reproducible results because lavaan introduces randomness into the calculations
LMed_fit <- lavaan::sem(LMedModel, data = d, se = "bootstrap", missing = "fiml")
LMed_Sum <- lavaan::summary(LMed_fit, standardized = T, rsq = T, ci = TRUE)
LMed_ParEsts <- lavaan::parameterEstimates(LMed_fit, boot.ci.type = "bca.simple",
standardized = TRUE)
LMed_Sum
## lavaan 0.6.16 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 7
##
## Number of observations 251
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 1000
## Number of successful bootstrap draws 1000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## central_sensitization ~
## trm_symp (b) 0.479 0.037 12.851 0.000 0.409 0.557
## strssf__ (c_p) 0.034 0.013 2.691 0.007 0.008 0.060
## trauma_symptoms ~
## strssf__ (a) 0.101 0.021 4.762 0.000 0.059 0.143
## Std.lv Std.all
##
## 0.479 0.645
## 0.034 0.134
##
## 0.101 0.296
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .trm_sym (tr_.) -0.179 0.362 -0.495 0.621 -0.866 0.539
## .cntrl_s (cn_.) 0.337 0.200 1.685 0.092 -0.059 0.715
## Std.lv Std.all
## -0.179 -0.184
## 0.337 0.468
##
## Variances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .centrl_snstztn 0.267 0.022 12.000 0.000 0.224 0.310
## .trauma_symptms 0.861 0.063 13.657 0.000 0.724 0.978
## Std.lv Std.all
## 0.267 0.514
## 0.861 0.913
##
## R-Square:
## Estimate
## centrl_snstztn 0.486
## trauma_symptms 0.087
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## indirect 0.048 0.011 4.248 0.000 0.027 0.073
## direct 0.034 0.013 2.689 0.007 0.008 0.060
## total_c 0.082 0.015 5.382 0.000 0.052 0.112
## Std.lv Std.all
## 0.048 0.191
## 0.034 0.134
## 0.082 0.325
LMed_ParEsts
## lhs op rhs label
## 1 central_sensitization ~ trauma_symptoms b
## 2 central_sensitization ~ stressful_life_events c_p
## 3 trauma_symptoms ~ stressful_life_events a
## 4 trauma_symptoms ~1 trauma_symptoms.mean
## 5 central_sensitization ~1 central_sensitization.mean
## 6 central_sensitization ~~ central_sensitization
## 7 trauma_symptoms ~~ trauma_symptoms
## 8 stressful_life_events ~~ stressful_life_events
## 9 stressful_life_events ~1
## 10 indirect := a*b indirect
## 11 direct := c_p direct
## 12 total_c := c_p+(a*b) total_c
## est se z pvalue ci.lower ci.upper std.lv std.all std.nox
## 1 0.479 0.037 12.851 0.000 0.407 0.555 0.479 0.645 0.645
## 2 0.034 0.013 2.691 0.007 0.009 0.060 0.034 0.134 0.047
## 3 0.101 0.021 4.762 0.000 0.058 0.142 0.101 0.296 0.104
## 4 -0.179 0.362 -0.495 0.621 -0.852 0.558 -0.179 -0.184 -0.184
## 5 0.337 0.200 1.685 0.092 -0.066 0.705 0.337 0.468 0.468
## 6 0.267 0.022 12.000 0.000 0.231 0.324 0.267 0.514 0.514
## 7 0.861 0.063 13.657 0.000 0.741 0.997 0.861 0.913 0.913
## 8 8.083 0.000 NA NA 8.083 8.083 8.083 1.000 8.083
## 9 16.717 0.000 NA NA 16.717 16.717 16.717 5.880 16.717
## 10 0.048 0.011 4.248 0.000 0.028 0.073 0.048 0.191 0.067
## 11 0.034 0.013 2.689 0.007 0.009 0.060 0.034 0.134 0.047
## 12 0.082 0.015 5.382 0.000 0.052 0.113 0.082 0.325 0.114
# Y = central sensitization, X = stressful life events, M = trauma symptoms, W = familismo
Combined <- '
#equations
trauma_symptoms ~ a1*stressful_life_events + a2*familism + a3*stressful_life_events:familism
central_sensitization ~ c_p1*stressful_life_events + c_p2*familism + c_p3*stressful_life_events:familism + b*trauma_symptoms
#intercepts
trauma_symptoms ~ trauma_symptoms.mean*1
central_sensitization ~ central_sensitization.mean*1
#means, variances of W for simple slopes
familism ~ familism.mean*1
familism ~~ familism.var*familism
#index of moderated mediation, there will be an a and b path in the product
#if the a and/or b path is moderated, select the label that represents the moderation
imm := a3*b
#Note that we first create the indirect product, then add to it the product of the imm and the W level
indirect.SDbelow := a1*b + imm*(familism.mean - sqrt(familism.var))
indirect.mean := a1*b + imm*(familism.mean)
indirect.SDabove := a1*b + imm*(familism.mean + sqrt(familism.var))
#direct effect is also moderated so calculate with c_p1 + c_p3
direct.SDbelow := c_p1 + c_p3*(familism.mean - sqrt(familism.var))
direct.Smean := c_p1 + c_p3*(familism.mean)
direct.SDabove := c_p1 + c_p3*(familism.mean + sqrt(familism.var))
'
set.seed(230925) #required for reproducible results because lavaan introduces randomness into the calculations
Combined_fit <- lavaan::sem(Combined, data = d, se = "bootstrap", missing = 'fiml', bootstrap = 1000)
## Warning in lav_partable_vnames(FLAT, "ov.x", warn = TRUE): lavaan WARNING:
## model syntax contains variance/covariance/intercept formulas
## involving (an) exogenous variable(s): [familism]; These variables
## will now be treated as random introducing additional free
## parameters. If you wish to treat those variables as fixed, remove
## these formulas from the model syntax. Otherwise, consider adding
## the fixed.x = FALSE option.
cFITsum <- lavaan::summary(Combined_fit, standardized = TRUE, rsq=T, ci=TRUE)
cParamEsts <- lavaan::parameterEstimates(Combined_fit, boot.ci.type = "bca.simple", standardized=TRUE)
cFITsum
## lavaan 0.6.16 ended normally after 17 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 13
##
## Number of observations 251
## Number of missing patterns 1
##
## Model Test User Model:
##
## Test statistic 946.622
## Degrees of freedom 2
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 1000
## Number of successful bootstrap draws 1000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## trauma_symptoms ~
## strss__ (a1) -0.031 0.089 -0.354 0.724 -0.205 0.150
## familsm (a2) -0.830 0.389 -2.132 0.033 -1.585 -0.082
## strs__: (a3) 0.029 0.023 1.292 0.196 -0.014 0.074
## central_sensitization ~
## strss__ (c_p1) 0.092 0.049 1.854 0.064 -0.006 0.193
## familsm (c_p2) 0.214 0.235 0.910 0.363 -0.254 0.694
## strs__: (c_p3) -0.016 0.013 -1.211 0.226 -0.043 0.010
## trm_sym (b) 0.461 0.042 10.986 0.000 0.382 0.544
## Std.lv Std.all
##
## -0.031 -0.070
## -0.830 -0.622
## 0.029 0.410
##
## 0.092 0.371
## 0.214 0.293
## -0.016 -0.413
## 0.461 0.843
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## .trm_sym (tr_.) 3.290 1.536 2.142 0.032 0.274 6.210
## .cntrl_s (cn_.) -0.395 0.886 -0.446 0.655 -2.205 1.364
## familsm (fml.) 3.673 0.060 60.861 0.000 3.555 3.785
## Std.lv Std.all
## 3.290 2.562
## -0.395 -0.563
## 3.673 3.820
##
## Variances:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## familsm (fml.) 0.924 0.081 11.469 0.000 0.769 1.075
## .trm_sym 0.763 0.060 12.821 0.000 0.628 0.863
## .cntrl_s 0.262 0.022 11.981 0.000 0.217 0.300
## Std.lv Std.all
## 0.924 1.000
## 0.763 0.463
## 0.262 0.532
##
## R-Square:
## Estimate
## trauma_symptms 0.537
## centrl_snstztn 0.468
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) ci.lower ci.upper
## imm 0.014 0.011 1.278 0.201 -0.007 0.035
## indirect.SDblw 0.022 0.015 1.494 0.135 -0.006 0.053
## indirect.mean 0.035 0.010 3.423 0.001 0.016 0.056
## indirect.SDabv 0.048 0.014 3.452 0.001 0.022 0.077
## direct.SDbelow 0.048 0.017 2.807 0.005 0.016 0.084
## direct.Smean 0.032 0.013 2.556 0.011 0.007 0.058
## direct.SDabove 0.016 0.019 0.873 0.383 -0.021 0.053
## Std.lv Std.all
## 0.014 0.346
## 0.022 0.917
## 0.035 1.263
## 0.048 1.608
## 0.048 -0.792
## 0.032 -1.205
## 0.016 -1.618
cParamEsts
## lhs op
## 1 trauma_symptoms ~
## 2 trauma_symptoms ~
## 3 trauma_symptoms ~
## 4 central_sensitization ~
## 5 central_sensitization ~
## 6 central_sensitization ~
## 7 central_sensitization ~
## 8 trauma_symptoms ~1
## 9 central_sensitization ~1
## 10 familism ~1
## 11 familism ~~
## 12 trauma_symptoms ~~
## 13 central_sensitization ~~
## 14 stressful_life_events ~~
## 15 stressful_life_events ~~
## 16 stressful_life_events:familism ~~
## 17 stressful_life_events ~1
## 18 stressful_life_events:familism ~1
## 19 imm :=
## 20 indirect.SDbelow :=
## 21 indirect.mean :=
## 22 indirect.SDabove :=
## 23 direct.SDbelow :=
## 24 direct.Smean :=
## 25 direct.SDabove :=
## rhs label
## 1 stressful_life_events a1
## 2 familism a2
## 3 stressful_life_events:familism a3
## 4 stressful_life_events c_p1
## 5 familism c_p2
## 6 stressful_life_events:familism c_p3
## 7 trauma_symptoms b
## 8 trauma_symptoms.mean
## 9 central_sensitization.mean
## 10 familism.mean
## 11 familism familism.var
## 12 trauma_symptoms
## 13 central_sensitization
## 14 stressful_life_events
## 15 stressful_life_events:familism
## 16 stressful_life_events:familism
## 17
## 18
## 19 a3*b imm
## 20 a1*b+imm*(familism.mean-sqrt(familism.var)) indirect.SDbelow
## 21 a1*b+imm*(familism.mean) indirect.mean
## 22 a1*b+imm*(familism.mean+sqrt(familism.var)) indirect.SDabove
## 23 c_p1+c_p3*(familism.mean-sqrt(familism.var)) direct.SDbelow
## 24 c_p1+c_p3*(familism.mean) direct.Smean
## 25 c_p1+c_p3*(familism.mean+sqrt(familism.var)) direct.SDabove
## est se z pvalue ci.lower ci.upper std.lv std.all std.nox
## 1 -0.031 0.089 -0.354 0.724 -0.210 0.137 -0.031 -0.070 -0.024
## 2 -0.830 0.389 -2.132 0.033 -1.589 -0.083 -0.830 -0.622 -1.839
## 3 0.029 0.023 1.292 0.196 -0.013 0.075 0.029 0.410 0.023
## 4 0.092 0.049 1.854 0.064 -0.006 0.193 0.092 0.371 0.131
## 5 0.214 0.235 0.910 0.363 -0.259 0.686 0.214 0.293 0.865
## 6 -0.016 0.013 -1.211 0.226 -0.043 0.010 -0.016 -0.413 -0.023
## 7 0.461 0.042 10.986 0.000 0.384 0.546 0.461 0.843 0.843
## 8 3.290 1.536 2.142 0.032 0.275 6.211 3.290 2.562 2.562
## 9 -0.395 0.886 -0.446 0.655 -2.222 1.342 -0.395 -0.563 -0.563
## 10 3.673 0.060 60.861 0.000 3.553 3.784 3.673 3.820 3.820
## 11 0.924 0.081 11.469 0.000 0.785 1.101 0.924 1.000 1.000
## 12 0.763 0.060 12.821 0.000 0.656 0.889 0.763 0.463 0.463
## 13 0.262 0.022 11.981 0.000 0.227 0.321 0.262 0.532 0.532
## 14 8.083 0.000 NA NA 8.083 8.083 8.083 1.000 8.083
## 15 20.272 0.000 NA NA 20.272 20.272 20.272 0.399 20.272
## 16 319.186 0.000 NA NA 319.186 319.186 319.186 1.000 319.186
## 17 16.717 0.000 NA NA 16.717 16.717 16.717 5.880 16.717
## 18 60.825 0.000 NA NA 60.825 60.825 60.825 3.405 60.825
## 19 0.014 0.011 1.278 0.201 -0.006 0.035 0.014 0.346 0.019
## 20 0.022 0.015 1.494 0.135 -0.007 0.052 0.022 0.917 0.034
## 21 0.035 0.010 3.423 0.001 0.016 0.056 0.035 1.263 0.053
## 22 0.048 0.014 3.452 0.001 0.022 0.077 0.048 1.608 0.073
## 23 0.048 0.017 2.807 0.005 0.016 0.084 0.048 -0.792 0.065
## 24 0.032 0.013 2.556 0.011 0.007 0.058 0.032 -1.205 0.042
## 25 0.016 0.019 0.873 0.383 -0.019 0.053 0.016 -1.618 0.019