g_REML.RdEstimates a design-comparable standardized mean difference effect size based on data from a multiple baseline design, using adjusted REML method as described in Pustejovsky, Hedges, & Shadish (2014). Note that the data must contain one row per measurement occasion per case.
g_REML(
m_fit,
p_const,
r_const,
X_design = model.matrix(m_fit, data = m_fit$data),
Z_design = model.matrix(m_fit$modelStruct$reStruct, data = m_fit$data),
block = nlme::getGroups(m_fit),
times = attr(m_fit$modelStruct$corStruct, "covariate"),
returnModel = TRUE
)Fitted model of class lme, with AR(1) correlation structure at level 1.
Vector of constants for calculating numerator of effect size.
Must be the same length as fixed effects in m_fit.
Vector of constants for calculating denominator of effect size.
Must be the same length as the number of variance component parameters in m_fit.
(Optional) Design matrix for fixed effects. Will be extracted from m_fit if not specified.
(Optional) Design matrix for random effects. Will be extracted from m_fit if not specified.
(Optional) Factor variable describing the blocking structure. Will be extracted from m_fit if not specified.
(Optional) list of times used to describe AR(1) structure. Will be extracted from m_fit if not specified.
(Optional) If true, the fitted input model is included in the return.
A list with the following components
p_beta | Numerator of effect size |
r_theta | Squared denominator of effect size |
delta_AB | Unadjusted (REML) effect size estimate |
nu | Estimated denominator degrees of freedom |
kappa | Scaled standard error of numerator |
g_AB | Corrected effect size estimate |
V_g_AB | Approximate variance estimate |
cnvg_warn | Indicator that model did not converge |
sigma_sq | Estimated level-1 variance |
phi | Estimated autocorrelation |
Tau | Vector of level-2 variance components |
I_E_inv | Expected information matrix |
Pustejovsky, J. E., Hedges, L. V., & Shadish, W. R. (2014). Design-comparable effect sizes in multiple baseline designs: A general modeling framework. Journal of Educational and Behavioral Statistics, 39(4), 211-227. doi:10.3102/1076998614547577
data(Laski)
Laski_RML <- lme(fixed = outcome ~ treatment,
random = ~ 1 | case,
correlation = corAR1(0, ~ time | case),
data = Laski)
summary(Laski_RML)
#> Linear mixed-effects model fit by REML
#> Data: Laski
#> AIC BIC logLik
#> 1048.285 1062.466 -519.1424
#>
#> Random effects:
#> Formula: ~1 | case
#> (Intercept) Residual
#> StdDev: 15.68278 13.8842
#>
#> Correlation Structure: AR(1)
#> Formula: ~time | case
#> Parameter estimate(s):
#> Phi
#> 0.252769
#> Fixed effects: outcome ~ treatment
#> Value Std.Error DF t-value p-value
#> (Intercept) 39.07612 5.989138 119 6.524498 0
#> treatmenttreatment 30.68366 2.995972 119 10.241637 0
#> Correlation:
#> (Intr)
#> treatmenttreatment -0.272
#>
#> Standardized Within-Group Residuals:
#> Min Q1 Med Q3 Max
#> -2.72642154 -0.69387388 0.01454473 0.69861200 2.14528141
#>
#> Number of Observations: 128
#> Number of Groups: 8
g_REML(Laski_RML, p_const = c(0,1), r_const = c(1,0,1), returnModel=FALSE)
#> Warning: 'g_REML()' is deprecated and may be removed in a later version of the package. Please use 'g_mlm()' instead.
#> est se
#> unadjusted effect size 1.465 0.299
#> adjusted effect size 1.405 0.286
#> degree of freedom 18.552
data(Schutte)
Schutte$trt.week <- with(Schutte, unlist(tapply((treatment=="treatment") * week,
list(treatment,case), function(x) x - min(x))) + (treatment=="treatment"))
Schutte$week <- Schutte$week - 9
Schutte_RML <- lme(fixed = fatigue ~ week + treatment + trt.week,
random = ~ week | case,
correlation = corAR1(0, ~ week | case),
data = subset(Schutte, case != 4))
summary(Schutte_RML)
#> Linear mixed-effects model fit by REML
#> Data: subset(Schutte, case != 4)
#> AIC BIC logLik
#> 875.9526 901.8978 -428.9763
#>
#> Random effects:
#> Formula: ~week | case
#> Structure: General positive-definite, Log-Cholesky parametrization
#> StdDev Corr
#> (Intercept) 9.783287 (Intr)
#> week 1.412030 0.811
#> Residual 5.421656
#>
#> Correlation Structure: AR(1)
#> Formula: ~week | case
#> Parameter estimate(s):
#> Phi
#> 0.3985953
#> Fixed effects: fatigue ~ week + treatment + trt.week
#> Value Std.Error DF t-value p-value
#> (Intercept) 50.29150 4.073625 121 12.345637 0.0000
#> week 0.20271 0.616194 121 0.328968 0.7427
#> treatmenttreatment -0.54218 1.751668 121 -0.309520 0.7575
#> trt.week -1.63198 0.655256 121 -2.490607 0.0141
#> Correlation:
#> (Intr) week trtmnt
#> week 0.883
#> treatmenttreatment -0.258 -0.200
#> trt.week -0.567 -0.596 -0.178
#>
#> Standardized Within-Group Residuals:
#> Min Q1 Med Q3 Max
#> -2.33194245 -0.46272153 -0.05348694 0.32878028 4.08749914
#>
#> Number of Observations: 136
#> Number of Groups: 12
Schutte_g <- g_REML(Schutte_RML, p_const = c(0,0,1,7), r_const = c(1,0,1,0,0))
#> Warning: 'g_REML()' is deprecated and may be removed in a later version of the package. Please use 'g_mlm()' instead.
summary(Schutte_g)
#> est se
#> sigma_sq 29.394 6.670
#> phi 0.399 0.132
#> case.var((Intercept)) 95.713 46.939
#> case.cov(week,(Intercept)) 11.209 6.372
#> case.var(week) 1.994 1.076
#> total variance 125.107 46.784
#> (Intercept) 50.291 4.074
#> week 0.203 0.616
#> treatmenttreatment -0.542 1.752
#> trt.week -1.632 0.655
#> treatment effect at a specified time -11.966 4.609
#> unadjusted effect size -1.070 0.495
#> adjusted effect size -1.013 0.469
#> degree of freedom 14.302
#> kappa 0.412
#> logLik -428.976