Changelog

Source: NEWS.md

JointAI (development version)

Small improvements

• Labelled variables (class "haven_labelled" from haven, or "labelled" from Hmisc and older versions of haven), e.g., from data imported from SPSS, can now be used (#17). Variables with value labels are converted to unordered factors (with a message); other labelled variables keep their values. Variable labels are kept. SPSS user-defined missing values are not treated as missing. Labelled variables in newdata of predict() are converted in the same way, using the categories of the fitted model.
• Clearer error messages: plot_imp_distr() applied to a fitted model now points to get_MIdat(); residuals() for JM and mlogit models says that residuals (not prediction) are not yet implemented; typo in the message for coxph models fixed.
• sum_duration(): the possible values of by are now listed in the function definition, with the new default "all" (NULL, the earlier default, still works); other values give an informative error instead of returning NULL.
• The deprecation warning of internal deprecated functions now names the function also when it is called via lapply() (it said “‘FUN’ is deprecated”), and the name of the replacement no longer has a trailing space.

Bug fixes

• A model type given in models for a covariate must now fit the level of the variable: with random effects (e.g., "lmm") for a variable that varies within the groups of a grouping level, without random effects (e.g., "lm") for a variable on the highest level. Before, a type without random effects for a level-one variable in a multilevel model silently dropped the intercept and the higher-level covariates from the linear predictor of that variable’s model; the other mismatches stopped with an error from JAGS. The new error message names the type that fits.
• Names of the regression coefficients in the output: summary() with a subset of the coefficients (subset = list(other = ...)) showed the values under the names of other coefficients (e.g. beta[3] as “(Intercept)”, with a warning “number of items to replace …”). traceplot(), densplot(), GR_crit() and MC_error() (models with one outcome) swapped names when the coefficients were not in the order of the coefficient table: with a subset in a different order, and by default in ordinal models with non-proportional effects and an interaction (e.g. o1 ~ c1 * C2 with nonprop = ~ c1 + C2, where c1:C2 was shown for a non-proportional effect of C2).
• list_models() describes the whole model (#113): it now lists the non-proportional effects of cumulative logit models, truncation, the random effects per grouping level with their covariance matrix and its prior (also for random effects in proportional hazards models), and the coefficients of multinomial mixed models per category (they were listed per design matrix, with the category names repeated). Small layout corrections.
• Models with random effects for a grouping variable within which the outcome does not vary (e.g., an outcome on the level of the subjects with (1 | id)) now give an informative error. They failed in JAGS, or, for multinomial mixed models, gave a model without covariates and random effects without any message.
• traceplot(), densplot() and predict(): sub-models can now be selected by name in outcome (e.g. outcome = "chol"), as documented; this failed with “‘start’ cannot be after ‘end’” (plots) or “EXPR must be a length 1 vector” (predict()). A value that is not a sub-model of the model now gives an informative error.
• monitor_params (in the *_imp() functions and add_samples()) combined with c(), e.g. c(analysis_main = TRUE, other = "tau_y"), now monitors the same parameters as with list(); the logical values were silently ignored, because c() turns them into character strings (#114). Unknown or unnamed elements of monitor_params give an informative error.
• subset in summary(), coef(), confint(), traceplot(), densplot(), GR_crit() and MC_error() combined with c() now selects the same parameters as with list(): e.g. c(analysis_main = FALSE, other = "beta[2]") included all parameters of the analysis model, because c() turns FALSE into "FALSE". Unknown or unnamed elements give an informative error.
• Misspelled arguments of the *_imp() functions (e.g. n.itre) are no longer ignored silently, but give a warning that names them (#114).
• JM_imp() with a survival part without covariates (e.g., Surv(time, status) ~ 1) gives an informative error; it stopped with “‘names’ attribute [2] must be the same length as the vector [0]”.
• Survival outcomes: Surv() with named arguments, e.g. Surv(time = futime, event = status) or Surv(futime, event = status), now gives the same model as Surv(futime, status); it failed with a syntax error in the JAGS model. Models with time-varying covariates no longer fail when the name of the survival time is part of the name of the id variable (e.g. t and patient).
• Informative errors instead of errors from JAGS: for a joint model with a time-varying unordered factor with more than two categories (its values between the observation times cannot be modelled yet), and for linear, generalized linear and multinomial (mixed) models without any regression coefficient (e.g., y ~ 0 + (1 | id)).
• survreg_imp(), and the survival part of JM_imp(), with more than one random effect on a level (e.g. (x | center)) failed in JAGS (“Unknown variable KinvD_…”), because the hyper-parameters of the prior of the random effects covariance matrix were missing from the data.
• auto_corr(), auto_corr_plot(), cross_corr(), cross_corr_plot(): the argument outcome now restricts the output to the selected sub-models, as was documented; it had only determined which coefficients were named by variable, while the others kept the name of their node (beta[3]). Its default changes from 1L to NULL, which includes all monitored parameters, as before. With parameters from several sub-models the coefficients are named "outcome: variable".
• coxph_imp() with time-varying covariates (timevar) and a list of formulas in which a survival formula comes after a formula for another outcome (e.g., list(Surv(...) ~ age, y ~ day, Surv(...) ~ x)) failed in JAGS, because the last observation carried forward was applied to the covariates of the first formulas instead of those of the survival formulas.
• Checks for an outcome that is used in more than one model formula:
  • Two outcomes of survreg_imp() with the same survival times (e.g., Surv(time, cause == 1) and Surv(time, cause == 2)) now give an error that says so, and that competing risks are not supported. (In coxph_imp() and JM_imp(), several survival outcomes can still share the survival time.)
  • An event indicator that is used in two survival outcomes, or a survival time or event indicator that is also the outcome of another model, now gives an error in all survival models.
  • A copy of an outcome (an outcome with identical values under another name) now gives an error that says so, instead of claiming that the variables were used in more than one formula. Identical values stored as integer and as double were not recognised.
  • Two outcomes of survreg_imp() whose event indicators have identical values (e.g., no censoring in either outcome) gave a false error.
• JM_imp(): exp() of a modelled time-varying covariate in the survival part (exp(x), or within I()) failed in JAGS with “Invalid vector argument to exp”, when using the association type "underl.value" (the default for, e.g., covariates with a Gaussian model).
• JM_imp(): specifications that are not the intended use of a joint model (#110). A longitudinal model whose outcome is not in the survival part failed with a cryptic error. It can now be fitted, with a warning that this outcome is assumed to be independent of the survival outcome (also for an outcome that is not used anywhere in a model that otherwise is a joint model). A model without a random effects structure was set up with one subject per row of the data, and time-varying covariates were used as baseline covariates, without error or warning: this is now an error, as is a model with neither a time-varying term in the survival part nor a longitudinal model, e.g., Surv(time, event) ~ age + (1 | id), which points to coxph_imp(). There are new warnings for a time-varying covariate in the survival part without a model in the formula (the automatically added model only has a random intercept), and for a model in which the only time-varying terms are effects of baseline covariates that change with time (e.g., age:day), which is not really a joint model.
• JM_imp(): when the survival part contained a function of several baseline variables (e.g., I(age * trig)), the association types were linked to the wrong variables. summary() and list_models() showed the wrong variable or NA, and an assoc_type given by the user was ignored or applied to another time-varying covariate, without error or warning. Results of such models with assoc_type fitted with earlier versions differ.
• JM_imp(): a function of several variables of which at least one is time-varying (e.g., I(age * chol)) in the survival part now gives an error that names the term and points to the interaction (age:chol). It stopped with “values must be length 1”. Such a model still cannot be fitted.
• JM_imp(): names in assoc_type that are not a time-varying term in the survival part now give an error that lists the possible names. They were silently ignored, e.g., c(chol = "obs.value") for a model with log(chol), where c("log(chol)" = "obs.value") is needed.
• survreg_imp(), coxph_imp() and JM_imp() now give an informative error for survival outcomes that are not right-censored. Left-censored outcomes (Surv(time, event, type = "left")) were accepted and fitted as if they were right-censored, without any message. Start-stop and interval-censored outcomes (e.g., Surv(start, stop, event)) failed with an error about competing risks.
• Factors with contr.sum or contr.helmert (set via options("contrasts")):
  • Incomplete factors with contr.sum stopped with an error. They now work, using the category number (not the label) in the JAGS model.
  • Complete factors with contr.sum did not get the documented coding (for more than two categories, only the first contrast variable was correct). They do now; the reference category is coded as -1 in all contrast variables, and the coefficients are named after the category they belong to (e.g., fb, fc instead of f1, f2). Results of such models fitted with earlier versions differ.
  • In the development version, any factor with these contrasts stopped with an error about the names of dummy variables. Complete factors with contr.helmert work again, with coefficients f1, f2.
• In mixed models with random effects on two or more levels, the intercept and the covariates of a higher level were silently left out of the model if they could not be in the means of the random effects of that level. Their coefficients were only sampled from the prior, and the other coefficients were estimated without these terms. They are now part of the linear predictor. This affected, for example, (1 | id) + (0 + time | center) and (0 + time | id) + (0 + time | center), for all types of mixed models and for time-varying covariates in JM_imp(); mlogitmm_imp() also with (1 | id) + (1 | center); and models of an outcome without random effects on a level that exists because of another outcome, e.g., list(y ~ z + (1 | id) + (1 | center), x ~ z + (1 | id)). Results of such models fitted with earlier versions differ.
• In models with three or more levels, an incomplete interaction of covariates from different levels above the lowest (e.g., a subject-level and a center-level covariate) was calculated with the value of the wrong unit: the center was looked up with the subject number instead of the row of that subject in the data, so that subjects were combined with the value of another center. Results of such models fitted with earlier versions differ.
• Non-proportional effects (nonprop) of ordinal outcomes in mixed models and joint models. Results of such models fitted with earlier versions differ.
  • If all non-proportional effects were on a higher level than the outcome (e.g., a baseline covariate with a longitudinal outcome), they were silently left out of the model. Their coefficients were only sampled from the prior. Effects from the higher level of a model with random effects on two levels (e.g., (1 | id) + (1 | center)) were always left out.
  • If all were on the level of the outcome, clmm_imp() stopped with an error.
  • For an ordinal time-varying covariate in a joint model (JM_imp()), they were left out of the part of the model that calculates the covariate at the quadrature points.
• JM_imp() stopped with the JAGS error “Unknown variable M_idgk” (or similar) when a time-varying covariate had a baseline covariate in its model that is not part of the means of the random effects, e.g., without random intercept: albumin ~ age + day + (0 + day | id).
• With several ordinal outcomes whose names start in the same way (e.g., O1 and O11), summary() and print() stopped with an error, and coef(), confint() and parameters() were wrong without any message: the thresholds of O11 were missing or had wrong names. Specifying, e.g., subset = list(other = "gamma_O1[1]") also selected parameters of the other outcome. With outcomes such as y and y.1, summary() listed sigma_y.1 for y as well.
• predict() with type = "response" rounded the predicted values of Poisson models (glm_imp(), glmm_imp()) to integers, also in the returned sample. It now gives the expected counts. The rounded values were also used for the fitted values, so that residuals() were wrong as well (NaN or very large values where the fitted value was rounded to 0). Results of such models fitted with earlier versions differ.
• residuals() of models with more than one model formula used the wrong type of residuals: by default only the first outcome got working residuals, the second got Pearson and the third response residuals. Now all get working residuals. In a named type, outcomes that are not named also get working residuals, and a single named type only applies to the named outcome.
• When the analysis model had no intercept (e.g., lm_imp(y ~ 0 + x1 + x2)) and no random effects, the models for the incomplete covariates (the imputation models) had no intercept either, e.g., x2 ~ 0 + x1. They now have an intercept (or, for ordinal covariates, thresholds), as when the analysis model has an intercept. Results of such models fitted with earlier versions differ.
• With scaled covariates (the default), centring the covariates changed some models, because no intercept could take up the shift: models without intercept (e.g., lm_imp(y ~ 0 + x1 + x2), which then had an intercept -sum(beta * mean(x))), random slopes without random intercept (e.g., (0 + time | id)), and interactions with a random-slope variable whose other variable had no main effect (e.g., y ~ x + x:time + (time | id)). These variables are now scaled without being centred. Results of such models fitted with earlier versions differ (#66).
• With scaled covariates (the default), the thresholds of ordinal models (clm_imp(), clmm_imp(), and ordinal covariate models) and the coefficients of the baseline hazard of proportional hazards and joint models (coxph_imp(), JM_imp()) were reported for centred covariates instead of on the scale of the original data, and the intercept of survreg_imp() models was not transformed back at all. These parameters are now transformed back like the intercept of other models (#62). Results of such models fitted with earlier versions differ in these parameters. Predictions of ordinal and proportional hazards models do not change; predictions of survreg_imp() models change (they used the untransformed intercept).
• The variance of a single scaled random effect (e.g., a random slope without random intercept, (0 + time | id)) had the same value in all iterations of a chain (that of the first iteration), so that GR_crit() was Inf (#67).
• extract_state() failed with an unclear error (“$ operator is invalid for atomic vectors”) for models that do not contain the JAGS model (fitted with n.adapt = 0). It now gives an informative error, as it does for objects that are not of class JointAI.
• With scaled covariates (the default), the coefficients of mixed models with an interaction between a covariate of a higher level (e.g., a baseline variable) and a variable with a random slope, e.g., lme_imp(y ~ trt * time + (time | id)), were transformed back incorrectly: the interaction, its main effects and the intercept (or the thresholds / baseline hazard) did not correspond to the fitted model. Predictions of such models change accordingly (#74).
• In mixed models with more than one grouping level, a fixed effect could be part of the mean of random effects on two levels, e.g., the effect of C1 in lme_imp(y ~ C1 * time + (time | id) + (C1 | center)) (mean of the random slope of C1 on the center level and of the random intercept on the id level), or the effect of time with random slopes of time on both levels. The effect was then counted twice, so that the estimated coefficient did not correspond to the model. It is now only part of the mean of the random effects of the highest level. Results of such models fitted with earlier versions differ (#76).
• With scaled covariates (the default), the random effects of mixed models with more than one grouping level were transformed back incorrectly when a covariate is in a cross-level interaction with a random-slope variable and its main effect is in the mean of random effects of a higher level, e.g., lme_imp(y ~ C1 * time + (time | id) + (C1 | center)) or lme_imp(y ~ Zc * time + (time | id) + (1 | center)) with Zc a center-level covariate: part of the effect of the covariate was in the random effects of the wrong level (or, in the first example, counted twice). This back-transformation needs the value of the covariate per group; if it has missing values, the random effects of these groups are only transformed back when the imputed values are monitored (monitor_params = c(imps = TRUE)), otherwise they are set to NA with a warning (#83).
• With scaled covariates (the default), the random effects of mixed models and their variance-covariance matrices were not always transformed back correctly: the random effects did not take into account the centring of covariates that is moved into the intercept (or thresholds / baseline hazard) (#70); for groups with fewer distinct observations than random effects, the random effects were wrong (#64); with rd_vcov = "full" and several outcomes, random effects and their covariance matrix were not transformed back at all (#63). Random effects and their covariance matrix are now also transformed back when a variable with a random slope has missing values; the warning that this was not possible is removed. When only some random effects or elements of their covariance matrix are monitored, the parameters needed for the back-transformation are added to the monitored parameters (only the other random effects of the selected groups, the covariance matrix, and the regression coefficients of the model), with a message (see monitor_params in ?model_imp).
• Mixed models with more than one grouping level and an interaction between a random-slope variable and a covariate of a higher grouping level, e.g., lme_imp(y ~ Zc * time + (time | id) + (1 | center)) with Zc a center-level covariate, could not be fitted (JAGS stopped with “Index out of range”). The covariate is now taken from the correct rows and is scaled (#77).
• With scaled covariates (the default), monitoring only some parameters failed with “subscript out of bounds” after sampling when the intercept of a covariate model (e.g., other = "alpha[1]"), the thresholds of an ordinal model (gamma_main, gamma_other) or the baseline hazard (basehaz) were monitored without the other regression coefficients of that model. These coefficients are now added to the monitored parameters (with a message), as was already the case for the main model(s). A monitored beta[...] is no longer monitored twice.
• The error message for variables of an unsupported class now names each variable correctly (it showed, e.g., “x1”, “x2” for a variable “x” with two classes) (#17).
• Model fitting failed with “invalid term in model formula” when the package logistf (or another package using formula.tools) was loaded (#15).
• Model fitting failed with “row names contain missing values” when variables had a “label” attribute (e.g., data imported from SPSS) or other attributes that are lost when the data are subset (#35).
• get_MIdat() failed with “object ‘.Random.seed’ not found” in an R session in which no random numbers had been drawn yet (e.g., a new session in plain R), also when seed was given. The same applied to the generation of initial values when inits is a function. If there was no .Random.seed before the call, none is left behind afterwards.
• Binomial mixed models with family = binomial(link = "cloglog") (and covariate models "glmm_binomial_cloglog") used a log link instead of the complementary log-log link. Results of such models fitted with earlier versions differ.
• Arguments that the *_imp() functions pass on to the underlying model function (e.g., trunc, custom) failed with “object not found” when their value was a local variable of a function from which *_imp() was called, or with “the … list contains fewer than … elements” when *_imp() was called through a function that passes on its ....
• predict() with newdata failed when newdata did not contain the outcome variable(s) (“replacement has 0 rows”), when a continuous covariate had only two distinct values in newdata (it was converted to a factor), and when a factor was given with only some of its levels, or as numeric or character. newdata is now converted according to the data used to fit the model. Proportional hazards models need the event time in newdata (also for type = "lp", the log hazard at that time); there is now an informative error if it is missing.
• scale_vars given as a character vector of variable names failed with “invalid argument type”, and scale_vars = TRUE switched scaling off. Now only the named variables are scaled, and TRUE is the same as the default.
• predict() for proportional hazards models (coxph_imp()) failed with “non-conformable arguments” when no variables were scaled (scale_vars = FALSE).
• Model fitting (and add_samples()) failed with “‘x’ must be an array of at least two dimensions” when the MCMC sample of a chain had only one iteration (e.g., n.iter = 1, or n.iter = 2 with thin = 2) and variables were scaled.
• Ordinal outcomes that were modelled with their own formula, but not as the first formula in the list (e.g., stage ~ age + (1 | id) in JM_imp(), or lm_imp(list(y ~ ..., O2 ~ ...))), had an intercept in addition to the thresholds, so that the model was not identified (#71). This was introduced in version 1.1.0.
• rd_vcov() failed with “subscript out of bounds” when the random-effects covariance matrices of covariate models were monitored (e.g., with monitor_params = c(other_models = TRUE)). It now returns them, together with those of the analysis model.
• MC_error(): the Monte Carlo error on the sampling scale (models fitted with keep_scaled_mcmc = TRUE) is now calculated for the selected nodes (subset) instead of all monitored nodes, with the same names as on the data scale; plot(MC_error(...), data_scale = FALSE) labelled the nodes wrongly or failed. Without the scaled sample, it now gives an informative error.
• print(summary()): the “Sample size per chain” was one too small when the difference between the last and the first iteration was a multiple of the thinning interval (e.g., 49 instead of 50 for iterations 110:600 with thinning interval 10).
• print(summary()): for thinned samples, the printed first and last iteration were the values of start and end as given, also when these were not stored iterations. They are now the first and last iteration used.
• A variable whose name was equal to the name of a dummy variable of a factor in the model (e.g., a variable B11 and the dummy of level 1 of a factor B1), or two factors with dummy variables of the same name, were merged into one variable, without error or warning. Model fitting now stops with an informative error (#75).
• Joint models (JM_imp()) with an interaction of a time-varying covariate in the survival part (e.g., age * chol, sex * chol or chol:day, with chol modelled over time) wrote invalid JAGS code (e.g., mu_age:chol), so that the model could not be fitted. Such an interaction is now the product of its elements: the time-varying covariate as specified by its association type, a baseline covariate with its value for the subject, and the time variable with its value at the event time (or the quadrature points). A covariate that is only part of an interaction now also gets an association type.
• Joint models (JM_imp()) with an interaction in the time-varying part of the survival model with an element that is not a column of the data used in the model stopped with uninformative errors (e.g., “I do not know the model type”NA”“) or, for an interaction of a baseline covariate with a function of time (e.g., age:ns(day, df = 3)), were fitted with values that did not vary over time. They now stop with an informative error. This concerns
  • an interaction with a function of the time variable without the main effect of that function (e.g., chol:ns(day, df = 3) without ns(day, df = 3)). As a workaround, the function can be added to the data with the argument auxvars (e.g., auxvars = ~ ns(day, df = 3)); this does not change the models for incomplete covariates.
  • an interaction of a factor with a variable without the main effect of that variable (e.g., sex + sex:day). R codes such an interaction with all categories of the factor (sexmale:day and sexfemale:day), which together contain a main effect of that variable. For the time variable (or a function of it), this main effect cannot (or hardly) be distinguished from the baseline hazard, and an interaction coded without the reference category (only sexfemale:day) is currently not possible.
• parameters() did not recognise regression coefficients that are monitored as single elements or ranges of a node (e.g., other = "beta[2]", or the coefficients that are added with scaled covariates when, e.g., a random effect is monitored): outcome and varname were NA and a range had one row. If no regression coefficients were monitored, the columns were in a different order.
• parameters() listed the regression coefficients of all models when all monitored parameters were regression coefficients (e.g., only beta monitored), also those of models that were not monitored.
• A joint model (JM_imp()) with a term in the survival part that involves nothing but the time variable (e.g., Surv(futime, status) ~ chol + day or Surv(futime, status) ~ chol + ns(day, df = 3) with timevar = "day") wrote invalid JAGS code. Since the effect of time is described by the baseline hazard, such a model now stops with an informative error.
• sum_duration(by = "run"): for models run in parallel, the duration per run was the longest chain over the adaptive and sampling phase together, instead of the sum of the longest chain per phase, so that the runs did not add up to the total duration.
• sum_duration() sorted the runs as text, so that with ten or more runs (e.g., after repeated add_samples()), “run 10” came before “run 2”. The runs are now kept in their numeric order.
• add_samples() failed with “attempt to set an attribute on NULL” for a model without MCMC sample (n.iter = 0) and the default add = TRUE. The new samples are now the samples of the object.
• add_samples() did not respect keep_scaled_mcmc = TRUE for a model without MCMC sample, and failed with “Start, end and thin incompatible with data” when adding to a scaled sample with thinning (thin > 1).
• add_samples() for a model without JAGS model (n.adapt = 0) now gives an informative error.
• An interaction with an element that is not a column of the data used in the model, because its main effect is not in the model (e.g., age:ns(day, df = 3) without ns(day, df = 3)), was placed on the level of the other element(s) (e.g., the level of age), also when it varies within that level. Its value was then that of the first observation of each group, without error or warning. Such an interaction is now placed on the level on which it varies.
• add_samples() kept the fitted values and residuals of the first run, also with add = FALSE, so that plot() used these old values and residuals() used the old fitted values. They are now calculated from the sample of the returned object (or removed if that is not possible).
• When add_samples() failed while re-creating the JAGS model (e.g., because of invalid parameter values), the console output stayed redirected to a temporary file, so that the error message and all later output were not shown.
• A joint model (JM_imp()) with a transformed time-varying covariate in the survival part (e.g., log(chol)) and the association type "obs.value" (specified via assoc_type = c("log(chol)" = "obs.value"), or the default, e.g., for a covariate with a Poisson model) applied the transformation twice, i.e., it used log(log(chol)), also in interactions. This could lead to an error in JAGS (“Invalid parent values”) or, without any error or warning (e.g., for I(chol^2)), to a wrong model. Results of such models fitted with earlier versions differ.
• predict() for models with more than one outcome:
  • Ordinal models used the thresholds of another ordinal variable when its name starts with the name of the outcome (e.g., O1 and O11), if these were monitored. The predicted probabilities were wrong, without an error, or the prediction failed; there were also additional categories named NA.
  • Proportional hazards models with several survival outcomes failed for type = "lp" and "risk" (“non-conformable arguments”), and the object had no fitted values. The coefficients of the baseline hazard of the other outcomes were used as well.
  • A single named type, e.g., type = c(b = "response"), applied to all outcomes and its name was not checked, and with a partly named vector the other outcomes got a recycled value. Now the named outcomes get the given type, the others "lp", and names that are not outcomes give an error.
• Fitted models with a time variable (timevar), auxiliary variables (auxvars), or a covariate that is only used inside a function (e.g., I(C2^2)) are much smaller when saved: they no longer carry a copy of the data and all design matrices, which was stored with an internal formula and did not show in object.size() (e.g., 4.7 MB instead of 0.7 MB for a joint model on PBC). Results are not affected.

JointAI 1.1.0

CRAN release: 2026-01-30

New features

• Added functions crosscorr(), crosscorr_plot(), autocorr() and autocorr_plot() to assess cross-correlation and autocorrelation of MCMC samples.

Small improvements

• all_vars(): update of the function; can now handle an unspecified number of input objects and can extract variable names from formulas and character strings that are valid variable names
• objects provided to the data argument are now converted with as.data.frame(), making it possible to provide a tibble
• Elements of a list supplied to the append_data_list argument will now overwrite existing elements of the data_list with the same name. This enables a wider range of changes to the automatically generated JAGS model.

Bug fixes

• Error when using the same function twice in a model formula with different number of variables (e.g., I(x^2) and I(a/b)) fixed.
• Fix in identifying the number of workers used in parallel computation with the future package.
• Fix summary() printout showing the wrong MCMC settings when using thinning.
• rd_vcov(): fixed issue resulting in error when extracting random effects variance-covariance matrix from fitted model in case of a multivariate mixed model.

---

JointAI 1.0.6

CRAN release: 2024-04-02

• clean-up of helper functions and additional unit tests
• fix typos in argument names in helpfiles
• fix documentation syntax (CRAN NOTEs)

---

JointAI 1.0.5

CRAN release: 2023-04-27

(update request by CRAN)

• use a wrapper for bs() to avoid warning in CMD check when run on older version of R where splines::bs() did not yet have the argument warn.outside.

---

JointAI 1.0.4

CRAN release: 2022-09-03

New features

• md_pattern() has an additional argument sort_columns to provide the option to switch off the sorting of columns by number of missing values.
• sum_duration(): new function to get sum of computational time across chains, phases, runs, …

Bug fixes

• formula() now also return the model formula when add_samples() is used.
• JM_imp() with ordinal longitudinal outcome is not using the correct version of the linear predictor in the quadrature procedure approximating the integral in the survival.
• Bug causing the wrong elements of the data matrix to be monitored when monitor_params(imps = TRUE) in survival models fixed.
• predDF(): bugfix for models including auxiliary variables (which were previously not included into the data)
  
• nonprop: bugfix for non-proportional effects in covariate models

Small improvements

• Use of add_samples() will now result in the call element of a JointAI object being a list and no longer a nested list.
• Re-structuring of the JAGS model code for proportional hazards models to reduce computational time.
• change in how parallel computation is done: does not rely on doFuture and foreach any more, only package future is required.
• comp_info element of fitted JointAI object has changed due to the changes in parallel computing, and computational time is not reported separately for the adaptive phase and the sampling phase, and separately per chain when parallel computation was used.
• predict() is now a lot faster for proportional hazards models.

---

JointAI 1.0.3

CRAN release: 2021-11-28

New features

• custom: new argument in the main analysis functions that allows the user to replace the JAGS syntax for sub-models with custom syntax. The argument expects a named list where the names are the names of the variables for which the custom model will be used. This feature requires some knowledge of JAGS syntax.
• append_data_list: new argument in the main analysis functions that allows the user to add elements to the list of data that is passed to JAGS. This may be necessary for custom sub-models.
• rd_vcov: new argument in the main analysis functions that allows the specification of the structure of the random effects variance-covariance matrices in (multivariate) mixed models (and joint models).

Minor improvements and bug fixes

• Bug fix in re-scaling of random effects variance covariance matrix in multi-level models with >2 levels where some levels have only random intercept.
• Bug fix in generating the names of random effect nodes to monitor when there are multiple analysis models and some do not have random effects.
• Bug fix in plot_all() not displaying the variable name in the message shown for character string variables.
• Bug fix in set_refcat(): wasn’t displaying the factor labels correctly
• print.JointAI: bugfix for multivariate mixed models with full variance-covariance matrix of the random effects that had the same number of random effects for each outcome, for which the printing of this matrix resulted in an error.
• Bug fix: when using a function of a variable as auxiliary this is now (again) correctly used as covariate in the linear predictor of covariate models (bug was introduced in commit 15014dcd)

---

JointAI 1.0.2

CRAN release: 2021-01-13

New features

• rd_vcov(): new function added to extract the random effect variance- covariance matrices (posterior means)

Minor improvements and bug fixes

• change in the way the model formula is contained in the model call, which should make it possible to call *_imp() functions from within another function that has the model formula as an argument.
• fixed issue that resulted in an error when the data was previously attached to the search path

---

JointAI 1.0.1

CRAN release: 2020-11-16

Minor improvements and bug fixes

• data_list: omit data matrix M_* from data_list if ncol == 0
• data_list: syntax to checking which pos_* to include can handle the case with multiple grouping variables being on the same, lowest level; before pos_* was excluded for only one of them
• random effects: it is now possible to use different grouping levels in different sub-models (when providing a list of model formulas)
• add_samples(): remove unnecessary call to doFuture::registerDoFuture()
• predDF() bug fix: the parameter for all methods is now called object
• list_models() now also works for errored JointAI objects
• Bug fix for models with an interaction between a repeatedly measured variable and a random slope variable. (The term was wrongly written in the mean structure for the random effect instead of the main linear predictor.)

---

JointAI 1.0.0

CRAN release: 2020-08-31

This version of JointAI contains some major changes. To extend the package it was necessary to change the internal structure and it was not possible to assure backward compatibility.

New features

New analysis model types

• betareg_imp(): beta regression
• betamm_imp(): beta mixed model
• lognorm_imp(): log-normal regression model
• lognormmm_imp(): log-normal mixed model
• mlogit_imp(): multinomial logit model
• mlogitmm_imp(): multinomial mixed model
• JM_imp(): joint model for longitudinal and survival data

Hierarchical models with multiple levels of grouping

It is now possible to fit hierarchical models with more than one level of grouping, with nested as well as crossed random effects (check the help file) of the main model function for details on how to specify such random effects structures.

This does also apply to survival models, i.e., it is possible to specify a random effects structure to model survival outcomes in data with a hierarchical structure, e.g., in a multi-centre setting.

Proportional hazards model with time-dependent covariates

coxph_imp() can now handle time-dependent covariates using last-observation-carried-forward. This requires to add (1 | <id variable>) to the model formula to identify which rows belong to the same subject, and to specify the argument timevar to identify the variable that contains the observation time of the longitudinal measurements.

Multivariate models

By providing a list of model formulas it is possible to fit multiple analysis models (of different types) simultaneously. The models can share covariates and it is possible to have the response of one model as covariate in another model (in a sequential manner, however, not circular).

Partial proportional odds models for ordinal responses

As before, proportional odds are assumed by default for all covariates of a cumulative logit model. The argument nonprop accepts a one-sided formula or a named list of one-sided formulas in which the covariates are specified for which non-proportional odds should be assumed.

Additionally, the argument rev is available to specify a vector of names of ordinal responses for which the odds should be inverted. For details, see the the help file.

Other new features

• The functions lmer_imp() and glmer_imp() are aliases for lme_imp() and glme_imp().
• All mixed models can be specified using the nlme type specification (using fixed and random) as well as using specification as in lme4 (with a combined model formula for fixed and random effects).
• The argument df_basehaz in coxph_imp() and JM_imp() allows setting the number of degrees of freedom in the B-spline basis used to model the baseline hazard.
• The global setting options("contrasts") is no longer ignored. For completely observed covariates, any of the contrasts available in base R are possible and options("contrasts") is used to determine which contrasts to use. For covariates with missing values, effect coding (contr.sum) and dummy coding (contr.treatment) are possible. This means that for completely observed ordered factors the default is contr.poly, but for incomplete ordered factors the default is contr.treatment.

Other changes

• Within summary(), the argument multivariate to the function GR_crit() is now set to FALSE to avoid issues. The multivariate version can still be obtained by using GR_crit() directly.
• The parameters of the baseline hazard for coxph_imp() and JM_imp() are monitored automatically when “analysis_main = TRUE”.
• The function get_models() is no longer exported because it now requires more input variables and is no longer convenient to use. To obtain the default specification of the model types use one of the main functions (*_imp()), set n.adapt = 0 (and n.iter = 0), and obtain the model types as <mymodel>$models.
• list_models() now returns information on all sub-models, including the main analysis models (previously it included only information on covariates).
• The output of parameters() has changed. The function returns a list of matrices, one per analysis model, with information on the response variable, response category, name of the regression coefficient and its associated covariate.
• The argument missinfo = TRUE in summar() adds information on the missing values in the data involved in a JointAI model (number and proportion of complete cases, number and proportion of missing values per variable).
• The argument ridge has changed to shrinkage. If shrinkage = "ridge", a ridge prior is used for all regression coefficients in all sub-models. If a vector of response variable names is provided to shinkage, ridge priors are only used for the coefficients in these models.
• default_hyperpars(): The default hyper-parameters for random effects are no longer provided as a function but more consistently with the hyper-parameters for other model parts, as a vector (within the list of all hyper-parameters).
• default_hyperpars(): the default number of degrees of freedom in the Wishart distribution used for the inverse of the random effects covariance matrix is now the number of random effects + 1 (was the number of random effects before).
• The effect of using a seed in JointAI is now only local, i.e., in functions in which set.seed() is called, the random number generator state .Random.seed before setting the seed is recorded, and re-set to that value on exiting the function.
• In predDF() the argument var has changed to vars and expects a formula. This extends the functionality of predDF() to let multiple variables vary.
• Some of the elements of a JointAI object have changed. Potentially relevant for users:
  • JointAI objects now contain information on the computational time and environment they were run in: the element comp_info contains the time-stamp the model was fitted, the duration of the computation, the JointAI version number and the R session info.
  • The JAGS model is stored as character string in the element jagsmodel.
• The arguments warn and mess now also affect the output of rjags.
• The doFuture package is used for parallel computation instead of doParallel. Parallel computation is specified by setting a future::plan() for how the “future” should be handled. As a consequence, the arguments parallel and n.cores are no longer used. Information on the setting that was used with regards to parallel computation is returned in a JointAI object via comp_list$future.

Bug fixes

• Scaling of continuous covariates is no longer done in the data, but in the JAGS model. This fixes the issue that previously the estimated variance parameters of continuous covariates and variance parameters of their random effects were incorrect (this issue only affected the covariate models, not the main analysis model).

---

JointAI 0.6.1

CRAN release: 2020-02-12

Bug fixes

• Fixed bug that messed up coefficients in clmm covariate model when there are no baseline covariates in the model
• enable newdata for predict() that does not contain the outcome variable
• add_samples(): calculation of end() of MCMC samples fixed

---

JointAI 0.6.0

CRAN release: 2019-08-31

Bug fixes

• bug in add_samples() when used in parallel with thinning fixed
• bug fixed that occurred when a complete longitudinal categorical variable was used in a model that did not contain any incomplete baseline variables
• bug-fix for monitoring random effects
• fixed typo in selecting parameters in Gamma models
• predict() can now handle newdata with missing outcome values; predicted values for cases with missing covariates are NA (prediction with incomplete covariates is planned to be implemented in the future)
• bug-fix for get_MIdat() and plot_imp_distr() when only one variable has missing values
• bug-fix for longitudinal model with interaction with random slope variable
• bug-fix for model with multiple longitudinal ordinal incomplete covariates (fixed wrong selection of columns of the design matrix of longitudinal covariates in these models)

Minor changes

• moved message about bug reports to startup
• enabled inverse link by adding max restriction
• “:” in factor labels are automatically replaced by “_”
• argument ncores has changed to n.cores for consistency with n.iter, n.chains, etc.
• coxph_imp() does no longer use a counting process implementation but uses the likelihood in JAGS directly via the zeros trick

New Features / Extensions

• predict() now has an argument length to change number of evaluation points
• summary(), predict(), traceplot(), densplot(), GR_crit(), MC_error() now have an argument exclude_chains that allows to specify chains that should be omitted
• citation() now refers to a manuscript on arXiv
• glmm_lognorm available to impute level-1 covariates with a log-normal mixed model
• methods residuals() and plot() available for (some of the) main analysis types (details see documentation)
• argument models added to get_models() so that the user can specify to also include models for complete covariates (which are then positioned in the sequence of models according to the systematic used in JointAI). Specification of a model not needed for imputation prints a notification.
• JointAI objects (most types) now also include residuals and fitted values (so far, only using fixed effects)

---

JointAI 0.5.2

CRAN release: 2019-06-06

Bug fixes

• Error message in print.JointAI fixed

---

JointAI 0.5.1

CRAN release: 2019-05-07

Bug fixes

• bug in ordinal models with only completely observed variables fixed (all necessary data is not passed to JAGS)
• enable thinning when using parallel sampling
• matrix Xl is no longer included in data_list when it is not used in the model
• bug-fix in subset when specified as vector
• bug-fix in ridge regression (gave an error message)
• bug-fix in recognition of binary factors that are coded as numeric and have missing values
• bug-fix in summary: range of iterations is printed correctly now when argument end is used
• bug-fix: error that occurred in re-scaling when reference category was changed is solved
• bug-fix in survival models: coding of censoring variable fixed

Minor changes

• summary() calls GR_crit() with argument autoburnin = FALSE unless specified otherwise via ...
• when inits is specified as a function, the function is evaluated and the resulting list passed to JAGS (previously the function was passed to JAGS)
• the example data simong and simWide have changed (more variables, less subjects)
• added check if there are incomplete covariates before setting imp_pars = TRUE (when user specified via monitor_params or subset)
• in survreg_imp the sign of the regression coefficient is now opposite to match the one from survreg

---

JointAI 0.5.0

CRAN release: 2019-03-08

Important

• the argument meth has changed to models

Bug fixes

• add_samples(): bug that copied the last chain to all other chains fixed
• bug-fix for the order of columns in the matrix Xc, so that specification of functions of covariates in auxiliary variables works better
• adding vertical lines to a densplot() issue (all plots showed all lines) fixed
• nested functions involving powers made possible
• typo causing issue in poisson glm and glme removed

Minor changes

• plot_all(), densplot(), and traceplot() limit the number of plots on one page to 64 when rows and columns of the layout are not user specified (to avoid the ‘figure margins too large’ error)
• change in longDF example data: new version containing complete and incomplete categorical longitudinal variables (and variable names L1 and L2 changed to c1 and c2)
• Some minor changes in notes, warnings and error messages
• The function list_impmodels() changed to list_models() (but list_impmodels() is kept as an alias for now)
• improved handling of functional forms of covariates (also in longitudinal covariates and random effects)

New Features / Extensions

• clm_imp() and clmm_imp(): new functions for analysis of ordinal (mixed) models
• It is now possible to impute incomplete longitudinal covariates (continuous, binary and ordered factors).
• coxph_imp(): new function to fit Cox proportional hazards models with incomplete (baseline) covariates
• Argument no_model allows to specify names of completely observed variables for which no model should be specified (e.g., “time” in a mixed model)
• Shrinkage: argument ridge = TRUE allows to use shrinkage priors on the precision of the regression coefficients in the analysis model
• plot_all() can now handle variables from classes Date and POSIXt
• new argument parallel allows different MCMC chains to be sampled in parallel
• new argument ncores allows to specify the maximum number of cores to be used
• new argument seed added for reproducible results; also a sampler (.RNG.name) and seed value for the sampler (.RNG.seed) are set or added to user-provided initial values (necessary for parallel sampling and reproducibility of results)
• plot_imp_distr(): new function to plot distribution of observed and imputed values

---

JointAI 0.4.0

CRAN release: 2018-12-04

Bug fixes

• RinvD is no longer selected to be monitored in random intercept model (RinvD is not used in such a model)
• fixed various bugs for models in which only the intercept is used (no covariates)

Minor changes

• summary(): reduced default number of digits
• continuous variables with two distinct values are converted to factor
• argument meth now uses default values if only specified for subset of incomplete variables
• get_MIdat(): argument minspace added to ensure spacing of iterations selected as imputations
• densplot(): accepts additional options, e.g., lwd, col, …
• list_models() replaces the function list_impmodels() (which is now an alias)

Extensions

• coef() method added for JointAI object and summary.JointAI object
• confint() method added for JointAI object
• print() method added for JointAI object
• survreg_imp() added to perform analysis of parametric (Weibull) survival models
• glme_imp() added to perform generalized linear mixed modelling
• extended documentation; two new vignettes on MCMC parameters and functions for after the model is estimated; added messages about coding of ordinal variables

---

JointAI 0.3.4

Bug fixes

• traceplot(), densplot(): specification of nrow AND ncol possible; fixed bug when only nrow specified

JointAI 0.3.3

Bug fixes

• remove deprecated code specifying contrast.arg that now in some cases cause error
• fixed problem identifying non-linear functions in formula when the name of another variable contains the function name

# JointAI 0.3.2
# JointAI 0.3.1
## Bug fixes * plot_all() uses correct level-2 %NA in title * simWide: case with no observed bmi values removed * traceplot(), densplot(): ncol and nrow now work with use_ggplot = TRUE * traceplot(), densplot(): error in specification of nrow fixed * densplot(): use of color fixed * functions with argument subset now return random effects covariance matrix correctly * summary() displays output with row name when only one node is returned and fixed display of D matrix * GR_crit(): Literature reference corrected * predict(): prediction with varying factor fixed * no scaling for variables involved in a function to avoid problems with re-scaling
## Minor changes * plot_all() uses xpd = TRUE when printing text for character variables * list_impmodels() uses line break when output of predictor variables exceeds getOption("width") * summary() now displays tail-probabilities for off-diagonal elements of D * added option to show/hide constant effects of auxiliary variables in plots * predict(): now also returns newdata extended with prediction

JointAI 0.3.0

CRAN release: 2018-08-14

Bug fixes

• monitor_params is now checked to avoid problems when only part of the main parameters is selected
• categorical imputation models now use min-max trick to prevent probabilities outside [0, 1]
• initial value generation for logistic analysis model fixed
• bug-fix in re-ordering columns when a function is part of the linear predictor
• bug-fix in initial values for categorical covariates
• bug-fix in finding imputation method when function of variable is specified as auxiliary variable

Minor changes

• md.pattern() now uses ggplot, which scales better than the previous version
• lm_imp(), glm_imp() and lme_imp() now ask about overwriting a model file
• analysis_main = T stays selected when other parameters are followed as well
• get_MIdat(): argument include added to select if original data are included and id variable .id is added to the dataset
• subset argument uses same logit as monitor_params argument
• added switch to hide messages; distinction between messages and warnings
• lm_imp(), glm_imp() and lme_imp() now take argument trunc in order to truncate the distribution of incomplete variables
• summary() now omits auxiliary variables from the output
• imp_par_list is now returned from JointAI models
• cat_vars is no longer returned from lm_imp(), glm_imp() and lme_imp(), because it is contained in Mlist$refs

Extensions

• plot_all() function added
• densplot() and traceplot() optional with ggplot
• densplot() option to combine chains before plotting
• example datasets NHANES, simLong and simWide added
• list_impmodels to print information on the imputation models and hyper-parameters
• parameters() added to display the parameters to be/that were monitored
• set_refcat() added to guide specification of reference categories
• extension of possible functions of variables in model formula to (almost all) functions that are available in JAGS
• added vignettes Minimal Example, Visualizing Incomplete Data, Parameter Selection and Model Specification

---

JointAI 0.2.0

CRAN release: 2018-07-05

Bug fixes

• md_pattern(): does not generate duplicate plot any more
• corrected names of imputation methods in help file
• scaling when no continuous covariates are in the model or scaling is deselected fixed
• initial value specification for coefficient for auxiliary variables fixed
• get_MIdat(): imputed values are now filled in in the correct order
• get_MIdat(): variables imputed with lognorm are now included when extracting an imputed dataset
• get_MIdat(): imputed values of transformed variables are now included in imputed datasets
• problem with non valid names of factor labels fixed
• data matrix is now ordered according to order in user-specified meth argument

Minor changes

• md.pattern(): adaptation to new version of md.pattern() from the mice package
• internally change all NaN to NA
• allow for scaling of incomplete covariates with quadratic effects
• changed hyperparameter for precision in models with logit link from 4/9 to 0.001

Extensions

• gamma and beta imputation methods implemented