Changes in Version 0.7.3
  - Made the backward-compatibility checks of summary() output robust to rounding-boundary flips: summary.mlVAR rounds to 3 decimals, and ~5e-5 cross-platform optimizer jitter can flip the last printed digit when a value sits near a rounding boundary (observed on CRAN win-builder and Debian: 0.370 vs 0.369). The summary columns are now compared against the fit object itself (rounded identically) instead of 3-dp literals; a separate check pins the Beta standard errors against the 0.6.1 reference. No package code was changed

Changes in Version 0.7.2
  - Relaxed the tolerance of the backward-compatibility test checks that involve iterative model fits (lmer/lm) from 1e-6 to 1e-3: optimizer termination differs slightly across platforms/BLAS builds (a ~5e-5 relative difference was observed on CRAN's win-builder), which is numerical jitter rather than a behavioral change. Deterministic data-generation checks keep the 1e-6 tolerance. No package code was changed

Changes in Version 0.7.1
  - Fixed mlVARsample() erroring with "unused arguments (originalData = ...)". It re-fits the model by passing the stored input$ of the original fit back to mlVAR(), but that list records more than mlVAR()'s own arguments (originalData, scale_means/scale_sds and scaled are kept for predict/residuals). Only the arguments mlVAR() accepts are now passed on.
  - Fixed the example in ?mlVAR0, which simulated its data with mlVARsim0() - a function that no longer exists. It now uses mlVARsim() and the current plot API.
  - estimator = "JAGS" now gives its informative "not implemented yet" error message; previously this check was unreachable because match.arg() rejected "JAGS" first with a generic 'arg should be one of' error
  - Calling mlVAR() without 'idvar' now gives an informative error ("'idvar' must be supplied; multilevel estimation requires multiple subjects. For single-subject VAR, see the graphicalVAR package."). Previously a stopifnot(!missing(idvar)) produced a cryptic error and was followed by a dead branch that would have assigned all rows to a single implicit subject - a path that fails inside lmer anyway, since random effects require at least two groups
  - Fixed the rank-deficiency check silently REORDERING the retained variables: when linearly dependent variables were dropped, the survivors were kept in QR-pivot order rather than input order, changing the node order in all outputs. The retained variables now keep their original input order
  - For contemporaneous = "fixed"/"unique" with the lmer estimator, output$contemporaneous is now NULL instead of being overloaded with the estimated Theta covariance object (no second-step lmer models exist on that path; the field otherwise holds a list of lmer fits and nothing in the package read the overloaded value)
  - Fixed mlVARsim() storing the wrong quantity in the simulated model's mu$SD: it recorded the mu_SD argument range (e.g. c(1, 1)) rather than the realised per-node standard deviations. mu$SD now holds sqrt(diag(Omega[1:nNode, 1:nNode])), i.e. one realised SD per node from the (possibly shrunk) random-effect covariance. Also added a maxAbsValue argument (default 100, unchanged behaviour) controlling the outlier threshold that triggers shrink-and-retry, corrected a stale "> 10" comment (the code tested > 100), and removed an unreachable rmvnorm() branch inside the lag > 0 data-generation loop
  - Fixed a parallel-cluster leak: with nCores > 1, the PSOCK cluster created by mlVAR() (and by the vendored parSim()) was only stopped inside the two-step contemporaneous branch, so contemporaneous = "fixed"/"unique" runs and any error path left worker processes and their connections open. The cluster is now unconditionally stopped via on.exit()
  - nCores now uses the requested number of worker processes; previously one fewer worker than requested was spawned (nClust <- nCores - 1)
  - Fixed temporal = "fixed" producing a Beta$SD array with a spurious extra lag slice: the array was sized from unique(model$lag), which includes an NA level for between-subject predictor rows, so Beta$SD had one more lag slice than Beta$mean (e.g. a lag-1 fit gave dim(Beta$SD) = 3 3 2 vs dim(Beta$mean) = 3 3 1). Beta$SD is now sized from the NA-stripped lags, matching Beta$mean. summary.mlVAR also gained a defensive length check so any future temporal-array size mismatch errors clearly instead of silently recycling
  - Performance: the "Computing random effects" stage of the two-step contemporaneous estimation (contemporaneous = "orthogonal"/"correlated") recomputed all subjects' posthoc residual covariance matrices on every iteration of the per-subject loop, i.e. O(N^2) covariance computations for N subjects. The list is now computed once before the loop (O(N)) and is named by subject ID so the subject lookup is explicit rather than positional. Results are identical
  - Fixed estimator = "lm" crashing for non-numeric (e.g. character) subject IDs: per-subject residuals were collected in a matrix together with the id column, coercing all residuals to character and breaking the contemporaneous estimation (cov() error for contemporaneous = "unique"/"fixed", cryptic random-effects error for "correlated"/"orthogonal"). Residuals are now stored in a type-preserving data frame; results for numeric IDs are unchanged
  - Fixed the vendored parSim() exclude argument being inverted: filter(!!!exprs) KEPT the rows matching the exclude conditions instead of dropping them. Exclude expressions are now combined with OR and negated, so a condition/rep row is removed if it matches any exclude expression (as documented). An empty or missing exclude leaves all conditions in place
  - Fixed mlVARsample() validating nSample against the number of subjects in the original object BEFORE dropping improper subjects (non-stationary Beta matrices, outlying contemporaneous variance). The check now runs after filtering and compares against the surviving subject count, with an error message reporting how many subjects survived versus how many were requested
  - mlVARsample(): removed a leftover commandArgs() call that read the R session's command-line arguments inside the function (unused), and wrapped the missingness sample size (pMissing * total_obs) in round() so a non-integer size is no longer passed to sample()
  - plot.mlVAR() now validates the 'lag' argument for temporal networks and the 'subject' argument (for temporal and contemporaneous networks) up front, giving an informative error (e.g. "'lag' exceeds the number of lags in the model (1).", or a message reporting the number of subjects) instead of a bare out-of-bounds array/list subscript error
  - Documented in ?mlVAR that the reported fixed-effect p-values are Wald-type tests based on a normal approximation (2 * (1 - pnorm(abs(estimate / SE)))) and should be treated as heuristic, particularly with few subjects where a t-distribution would be more appropriate. No computation was changed
  - print() for mlVARcompare objects now handles ties in the "Best fitting model" AIC/BIC lines: when several models share the minimum criterion the line is suffixed with " (tie)". Output is otherwise unchanged
  - Internal: removed an unused windowSize variable in movingWindow() and clarified in a comment that maxEffects counts lagged predictors only (auto-regressions are always included and not counted), matching Stepwise(); behaviour is unchanged

Changes in Version 0.6.1
  - Fixed resimulate(): person-specific means are now computed as stationary means (I - B_between)^{-1} * alpha_i rather than using the raw lmer intercept, which was incorrect when between-person predictors are uncentered
  - Added resimulate() arguments: nTime for simulating custom-length time series per person, keep_missing to toggle whether the original variable-level missingness pattern is applied, and variance ("model" or "empirical") to choose between fully model-implied innovation covariance or model partial correlation structure scaled by empirical residual SDs
  - Added full_detrend argument to mlVAR() for removing systematic occasion effects (e.g., time-of-day trends) before estimation

Changes in Version 0.6
  - Added mlGGM() function for multi-level Gaussian Graphical Model estimation using single-step nodewise regression, simultaneously estimating within-cluster and between-cluster partial correlation networks from cross-sectional multilevel data
  - Added predict() S3 method for mlVAR objects: returns fitted values, residuals, and observed data aligned to original input data (with NAs for unpredictable rows). Supports scale_back argument to return values on original data scale, include_ids to include id/day/beep columns, and newdata for out-of-sample predictions. Works with both lmer and lm estimators.
  - Added resimulate() S3 method for mlVAR objects: generates simulated data from a fitted model using person-specific parameters for posterior predictive checks
  - Rewrote residuals() S3 method for mlVAR objects as a wrapper around predict(), now supporting scale_back and include_ids arguments
  - Added .groups = "drop" to all grouped summarise/summarize calls to suppress dplyr 1.0.0+ lifecycle warnings
  - Fixed bug in Stepwise(): selected wrong model at each iteration step
  - Fixed bug in importMplus: temporal fixed effects were only extracted for the first outcome variable
  - Fixed bug in Mplus correlation samples: all samples were computed from the 3rd sample instead of each respective sample
  - Fixed bug in mlVARsample: subjects with improper or non-stationary models were not correctly excluded from simulation
  - Added iteration cap (1000) to mlVARsim repeat loop to prevent infinite loops when stable parameters cannot be generated
  - Minor: avoided redundant eigendecomposition in forcePositive()

Changes in Version 0.5.5
  - Fixed R CMD check NOTE: no visible binding for global variable 'id' in movingWindow
  - Replaced deprecated dplyr functions: summarize_each/funs() with across(), summarise_each_() with across(), filter_() with filter()
  - Removed library() calls from parSim.R
  - Removed stringsAsFactors argument from cbind() calls in NodeWise.R (had no effect)
  - Replaced plyr::ddply and plyr::join with dplyr equivalents; removed plyr dependency
  - Fixed unreachable code after stop() in mlVAR (changed to warning)
  - Fixed duplicate return statement in randomEffects
  - Cleaned up NAMESPACE: removed duplicate export and unused plyr imports
  - Minimum dplyr version bumped to >= 1.0.0

Changes in Version 0.5.4
  - mlVAR with lmer estimation now returns 'step1_residuals' in the output

Changes in Version 0.5.3
  - Random effect SDs of Gamma_theta stored in $results$Gamma_Theta$SD are now based on analytric standard deviations

Changes in Version 0.5.2
  - Removed deprecated function mlVARsim0

Changes in Version 0.5.1
  - Fixed remaining deprecated dplyr functions

Changes in Version 0.5
  - The 'mlVARsample' function has been added to mlVAR
  - Added Myrthe Veenman to contributor list
  - Fixed a bug where contemporaneous standard deviations were reported as variances instead of standard deviations
  - Fixed a bug with the beepvar argument
  - Replaced deprecated dplyr functions
  - Added a warning for when a beep is used multiple times
  - The 'nonsig' argument in the plot method now defaults to 'show' when SD=TRUE
  - Fixed a bug in the summary method when fixed effects estimation was used

Changes in version 0.4.3
  o mlVAR now issues a warning when < 20 observations per subject are used
  o Fixed a bug with 'lmerResults2'
  o Now suppressing warnings and messages from lmer
  o Added a progress bar for computing random effects

Changes in version 0.4.2
  o Contemporaneous multi-level models are now returned in the output

Changes in version 0.4.1
  o mlVAR now uses correlations of residuals as estimate for the contemporaneous correlation matrix (not partial) if estimated inverse covariance matrix is not properly invetable
  o Added mlVARsample function to run a simulation study 
given a mlVAR object.
  o Fixed a bug with estimator = "mPlus"
  o mlVAR now gives a warning when between-subject networks could not be computed, rather than breaking with an uninformative error.

Changes in version 0.4
  o Added AR argument to mlVAR to fit AR models only
  o estimator = "Mplus" is now supported! Requires Mplus 8 to be installed.
    o Several arguments have been added to mlVAR to handle Mplus estimation

Changes in version 0.3.3
  o The plot method for mlVAR sim objects now uses nonsig = "show"
  o plot method now uses nonsig = "show" by default!
  o Summary method now shows p-values for contemporaneous effects
  o Several small bugfixes

Changes in version 0.3.1
  o The 'partial' argument in 'plot.mlVAR' now defaults to TRUE
  o Added 'contemporaneous' argument to mlVAR
  o Added 'lm' estimator for fitting unique VAR models per subject
  o Added 'rule' argument to plot.mlVAR to set the rule of choosing significance in nodewise GGM estimation

Changes in version 0.3
  o Complete rework of package! 
    o mlVAR, mlVARsim, and relevant methods have been completely rewritten
    o Now support contemporaneous effects and between-subjects effects
    o Old functions are now labeled mlVAR0, mlVARsim0, etcetera