Matching longitudinal methodology models with complex sampling design. It fits fixed and random effects models and covariance structured models so far. It also provides tools to perform statistical tests considering these specifications as described in : Pacheco, P. H. (2021). "Modeling complex longitudinal data in R: development of a statistical package." < https://repositorio.ufjf.br/jspui/bitstream/ufjf/13437/1/pedrohenriquedemesquitapacheco.pdf>.
Also know as: 'clm'(complex longitudinal model). Reference the scientific article: https://repositorio.ufjf.br/jspui/bitstream/ufjf/13437/1/pedrohenriquedemesquitapacheco.pdf
Matching longitudinal methodology models with complex sampling design. It fits fixed and random effects models and covariance structured models so far. It also provides tools to perform statistical tests considering these specifications.
The mmcsd package is in its initial development stage and does not yet have a version available on CRAN.
However, the development version of the package is now available for download and can be installed via the repository on GitHub with:
# install.packages("devtools")
devtools::install_github("phmpacheco-ufjf/mmcsd")
There are a few major functions available for usage. Here’s an example using BHPS survey data.
First, we fit the fixed and random effects model considering the survey sampling design. The sigma argument can be: identity, exchangeable, autorregressive or a custom square matrix.
library(mmcsd)
fit <- mmcsd(
score ~ wave + ageg + ecacg + qualifg,
waves = wave, ids = id,
weights = weight, stratum = strata, cluster = cluster,
data = example_data, sigma = "exchangeable"
)
summary(fit)
#> Call:
#> mmcsd(formula = score ~ wave + ageg + ecacg + qualifg, waves = wave,
#> ids = id, weights = weight, stratum = strata, cluster = cluster,
#> data = example_data, sigma = "exchangeable")
#>
#> Coefficients:
#>
#> Estimate Std. Error t value P(t>|t|)
#> ------------ --------- ----------- --------- --------- ----
#> (Intercept) 22.29 0.30 74.0021 0.0000 ***
#> wave -0.04 0.02 -2.7355 0.0062 **
#> ageg2 -0.65 0.21 -3.0434 0.0023 **
#> ageg3 -0.82 0.24 -3.4073 0.0007 ***
#> ageg4 -1.01 0.28 -3.6672 0.0002 ***
#> ecacg2 -0.87 0.11 -8.0011 0.0000 ***
#> ecacg3 -0.76 0.17 -4.4817 0.0000 ***
#> ecacg4 0.21 0.24 0.8678 0.3855
#> ecacg5 -2.03 0.14 -14.2889 0.0000 ***
#> qualifg2 -0.52 0.25 -2.0632 0.0391 *
#> qualifg3 -0.66 0.26 -2.5131 0.0120 *
#> qualifg4 -0.50 0.25 -2.0176 0.0436 *
#> qualifg5 -1.22 0.26 -4.6979 0.0000 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Number of individuals: 1340
#> Number of waves: 5
#>
#> Sampling Design:
#> Sample Weights: Yes
#> Stratification: Yes
#> Clusterization: Yes
#>
#> Covariance Matrix Provided:
#> Type: exchangeable
#> rho: 0.567
#>
#> T1 T2 T3 T4 T5
#> --- ------ ------ ------ ------ ------
#> T1 1.000 0.567 0.567 0.567 0.567
#> T2 0.567 1.000 0.567 0.567 0.567
#> T3 0.567 0.567 1.000 0.567 0.567
#> T4 0.567 0.567 0.567 1.000 0.567
#> T5 0.567 0.567 0.567 0.567 1.000
#> ---
#>
#> Sigma Estimated:
#>
#> T1 T3 T5 T7 T9
#> --- -------- -------- -------- -------- --------
#> T1 13.4627 7.6690 7.1577 6.3445 6.0808
#> T3 7.6690 11.5988 7.8561 6.5722 6.5063
#> T5 7.1577 7.8561 12.5644 7.9730 7.4928
#> T7 6.3445 6.5722 7.9730 11.9087 7.9353
#> T9 6.0808 6.5063 7.4928 7.9353 11.9837
#> ---
Then, we fit the covariance structured model. But first, we must choose the covariance structure we want to fit. Here’s an example for the Uniform Correlation Model - UCM:
sigmaThetaExpr_viewer("UCM", 5)
#>
#>
#> Table: sigma2u sigma2v
#>
#> | | T1 | T2 | T3 | T4 | T5 |
#> |:--|:-----------------:|:-----------------:|:-----------------:|:-----------------:|:-----------------:|
#> |T1 | sigma2u + sigma2v | sigma2u | sigma2u | sigma2u | sigma2u |
#> |T2 | sigma2u | sigma2u + sigma2v | sigma2u | sigma2u | sigma2u |
#> |T3 | sigma2u | sigma2u | sigma2u + sigma2v | sigma2u | sigma2u |
#> |T4 | sigma2u | sigma2u | sigma2u | sigma2u + sigma2v | sigma2u |
#> |T5 | sigma2u | sigma2u | sigma2u | sigma2u | sigma2u + sigma2v |
Finally, we fit the model using our previous random effects fitted model.
fitTheta <- cov_mmcsd(fit,
fittingType = "PML", sigmaThetaExpr = "UCM",
optimParams = list(par = c(7, 5))
)
summary(fitTheta)
#> Call:
#> cov_mmcsd(fit = fit, fittingType = "PML", sigmaThetaExpr = "UCM",
#> optimParams = list(par = c(7, 5)))
#>
#> Fitting Function Type: PML
#> Sigma Theta Type: UCM
#>
#> Covariance Params: sigma2u sigma2v
#>
#> Coefficients:
#>
#> Estimate Std. Error
#> -------- --------- -----------
#> sigma2u 7.1588 0.3880
#> sigma2v 5.1449 0.1532
#> ---
#>
#> Optimization Details:
#>
#> Converged: TRUE
#> Optim Method: L-BFGS-B
#> Iteractions Number - Function: 8
#> Iteractions Number - Gradient: 8
#> Optim Function Value: 0.1641
#>
#> Start Value Lower Upper
#> -------- ------------ ------ ------
#> sigma2u 7 - -
#> sigma2v 5 - -
#>
#>
#> Sigma Theta Matrix:
#>
#> -------- -------- -------- -------- --------
#> 12.3037 7.1588 7.1588 7.1588 7.1588
#> 7.1588 12.3037 7.1588 7.1588 7.1588
#> 7.1588 7.1588 12.3037 7.1588 7.1588
#> 7.1588 7.1588 7.1588 12.3037 7.1588
#> 7.1588 7.1588 7.1588 7.1588 12.3037
#> -------- -------- -------- -------- --------
#>
#>
#> Goodness of Fit Measures:
#>
#> RMR: 0.8796
#> AGFI: 0.9225