Model Specification

Interaction with multiple variables

Interactions between multiple variables can be specified using parentheses:

mod2a <- glm_imp(educ ~ gender * (age + smoke + creat),
                 data = NHANES, family = binomial(), n.adapt = 0)

The function parameters() returns a matrix off all parameters that are specified to be followed (column coef) and, for regression coefficients, the name of the variable the coefficient relates to (varname), the outcome variable of the respective model outcome. For multinomial models, which have multiple linear predictors, the column outcat identifies the category of the outcome the parameters refer to.

We use the function parameters() here and in other vignettes to demonstrate the effect that different model specifications have.

parameters(mod2a)
#> 
#> Note: "mod2a" does not contain MCMC samples.
#>    outcome outcat                   varname     coef
#> 1     educ   <NA>               (Intercept)  beta[1]
#> 2     educ   <NA>              genderfemale  beta[2]
#> 3     educ   <NA>                       age  beta[3]
#> 4     educ   <NA>               smokeformer  beta[4]
#> 5     educ   <NA>              smokecurrent  beta[5]
#> 6     educ   <NA>                     creat  beta[6]
#> 7     educ   <NA>          genderfemale:age  beta[7]
#> 8     educ   <NA>  genderfemale:smokeformer  beta[8]
#> 9     educ   <NA> genderfemale:smokecurrent  beta[9]
#> 10    educ   <NA>        genderfemale:creat beta[10]
#> 11   creat   <NA>               (Intercept) alpha[1]
#> 12   creat   <NA>              genderfemale alpha[2]
#> 13   creat   <NA>                       age alpha[3]
#> 14   creat   <NA>               smokeformer alpha[4]
#> 15   creat   <NA>              smokecurrent alpha[5]
#> 16   smoke   <NA>              genderfemale alpha[6]
#> 17   smoke   <NA>                       age alpha[7]

Higher level interactions

To specify interactions of a given level, ^ can be used:

# all two-way interactions:
mod2b <- glm_imp(educ ~ gender + (age + smoke + creat)^2,
                 data = NHANES, family = binomial(), n.adapt = 0)

parameters(mod2b)
#>    outcome outcat            varname     coef
#> 1     educ   <NA>        (Intercept)  beta[1]
#> 2     educ   <NA>       genderfemale  beta[2]
#> 3     educ   <NA>                age  beta[3]
#> 4     educ   <NA>        smokeformer  beta[4]
#> 5     educ   <NA>       smokecurrent  beta[5]
#> 6     educ   <NA>              creat  beta[6]
#> 7     educ   <NA>    age:smokeformer  beta[7]
#> 8     educ   <NA>   age:smokecurrent  beta[8]
#> 9     educ   <NA>          age:creat  beta[9]
#> 10    educ   <NA>  smokeformer:creat beta[10]
#> 11    educ   <NA> smokecurrent:creat beta[11]
#> 12   creat   <NA>        (Intercept) alpha[1]
#> 13   creat   <NA>       genderfemale alpha[2]
#> 14   creat   <NA>                age alpha[3]
#> 15   creat   <NA>        smokeformer alpha[4]
#> 16   creat   <NA>       smokecurrent alpha[5]
#> 17   smoke   <NA>       genderfemale alpha[6]
#> 18   smoke   <NA>                age alpha[7]

# all two- and three-way interactions:
mod2c <- glm_imp(educ ~ gender + (age + smoke + creat)^3,
                 data = NHANES, family = binomial(), n.adapt = 0)

parameters(mod2c)
#>    outcome outcat                varname     coef
#> 1     educ   <NA>            (Intercept)  beta[1]
#> 2     educ   <NA>           genderfemale  beta[2]
#> 3     educ   <NA>                    age  beta[3]
#> 4     educ   <NA>            smokeformer  beta[4]
#> 5     educ   <NA>           smokecurrent  beta[5]
#> 6     educ   <NA>                  creat  beta[6]
#> 7     educ   <NA>        age:smokeformer  beta[7]
#> 8     educ   <NA>       age:smokecurrent  beta[8]
#> 9     educ   <NA>              age:creat  beta[9]
#> 10    educ   <NA>      smokeformer:creat beta[10]
#> 11    educ   <NA>     smokecurrent:creat beta[11]
#> 12    educ   <NA>  age:smokeformer:creat beta[12]
#> 13    educ   <NA> age:smokecurrent:creat beta[13]
#> 14   creat   <NA>            (Intercept) alpha[1]
#> 15   creat   <NA>           genderfemale alpha[2]
#> 16   creat   <NA>                    age alpha[3]
#> 17   creat   <NA>            smokeformer alpha[4]
#> 18   creat   <NA>           smokecurrent alpha[5]
#> 19   smoke   <NA>           genderfemale alpha[6]
#> 20   smoke   <NA>                    age alpha[7]

In JointAI, interactions between any variables, observed or incomplete, variables on different levels of a hierarchical structure, can be handled. When an incomplete variable is involved, the interaction term is re-calculated within each iteration of the MCMC sampling, using the imputed values from the current iteration.

It is important not to pre-calculate interactions with incomplete variables as extra variables in the dataset, but to specify them in the model formula. Otherwise, imputed values of the main effect and interaction term will not match, and results may be incorrect.

What happens inside JointAI?

When a model formula includes a function of a complete or incomplete variable, the main effect of that variable is automatically added as an auxiliary variable. (For more info on auxiliary variables, see the section “Auxiliary variables”.) In the linear predictors of models for covariates, usually, only the main effects are used.

In mod3b from above, for example, the spline of age is used as predictor for SBP, but in the imputation model for bili, age enters with a linear effect.

list_models(mod3b, priors = FALSE, regcoef = FALSE, otherpars = FALSE)
#> Linear model for "SBP" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), ns(age, df = 2)1, ns(age, df = 2)2, genderfemale, I(bili^2), I(bili^3) 
#> 
#> 
#> Linear model for "bili" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale

The function list_models() prints information on all sub-models specified in a JointAI object. This includes the model(s) specified by the user via the model formula and all models for covariates that JointAI has specified automatically. Since here we are only interested in the predictor variables, we suppress printing of information on prior distributions, regression coefficients and other parameters by setting priors, regcoef and otherpars to FALSE.

When a function of a variable is specified as an auxiliary variable, this function is used (as well) in the models for covariates. For example, in mod3e, waist circumference (WC) is not part of the model for SBP, and the auxiliary variable I(WC^2) is used in the linear predictor of the imputation model for bili:

mod3e <- lm_imp(SBP ~ age + gender + bili, auxvars = ~ I(WC^2), data = NHANES)

list_models(mod3e, priors = FALSE, regcoef = FALSE, otherpars = FALSE)
#> Linear model for "SBP" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale, bili 
#> 
#> 
#> Linear model for "bili" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale, I(WC^2) 
#> 
#> 
#> Linear model for "WC" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale

When WC is used as predictor variable in the main model and a function of WC is specified as auxiliary variable, both variables are used as predictors in the imputation models:

mod3f <- lm_imp(SBP ~ age + gender + bili + WC, auxvars = ~ I(WC^2), data = NHANES)

list_models(mod3f, priors = FALSE, regcoef = FALSE, otherpars = FALSE)
#> Linear model for "SBP" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale, bili, WC 
#> 
#> 
#> Linear model for "bili" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale, WC, I(WC^2) 
#> 
#> 
#> Linear model for "WC" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale

When a function of a covariate is used in the linear predictor of the analysis model, and that function should also be used in the linear predictor of imputation models, it is necessary to also include that function in the argument auxvars:

mod3g <- lm_imp(SBP ~ age + gender + bili + I(WC^2), auxvars = ~ I(WC^2), data = NHANES)

list_models(mod3g, priors = FALSE, regcoef = FALSE, otherpars = FALSE)
#> Linear model for "SBP" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale, bili, I(WC^2) 
#> 
#> 
#> Linear model for "bili" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale, WC, I(WC^2) 
#> 
#> 
#> Linear model for "WC" 
#>    family: gaussian 
#>    link: identity 
#> * Predictor variables:
#>   (Intercept), age, genderfemale

Incomplete variables are always imputed on their original scale, i.e.,

• in mod3b the variable bili is imputed and the quadratic and cubic versions calculated from the imputed values.
• Likewise, creat and albu in mod3c are imputed separately, and I(creat/albu^2) calculated from the imputed (and observed) values.

Important:
When different transformations of the same incomplete variable are used in one model, it is strongly discouraged to calculate these transformations beforehand and to supply them as separate variables. The same is the case for interactions.
If, for example, a model formula contains both x and x2 (where x2 = x^2), they are treated as separate variables and imputed with different models. Imputed values of x2 are thus not equal to the square of imputed values of x. Instead, x + I(x^2) should be used in the model formula. Then, only x is imputed and used in the linear predictor of models for other incomplete variables, and x^2 is calculated from the imputed values of x.

Functions with restricted support

When a function has restricted support, e.g., log(x) is only defined for x > 0, the model used to impute x needs to comply with these restrictions. This can either be achieved by truncating the distribution assumed for x, using the argument trunc, or by specifying a model for x that meets the restrictions. For more information on imputation methods (models for covariates), see the section Imputation model types.

Example:
When using a log() transformation for the covariate bili, we can either use the default model for continuous variables, "lm", a linear model, i.e., assuming a normal distribution and truncate this distribution by specifying trunc = list(bili = c(<lower>, <upper>)) (where the lower and upper limits are the smallest and largest allowed values) or choose a model (using the argument models; more details see the section on covariate model types) that only imputes positive values such as a log-normal distribution ("glm_lognorm") or a Gamma distribution (e.g., "glm_gamma_log"):

# truncation of the distribution of  bili
mod4a <- lm_imp(SBP ~ age + gender + log(bili) + exp(creat),
                trunc = list(bili = c(1e-5, NA)), data = NHANES)

# log-normal model for bili
mod4b <- lm_imp(SBP ~ age + gender + log(bili) + exp(creat),
                models = c(bili = 'lognorm', creat = 'lm'), data = NHANES)

# gamma model with log-link for bili
mod4c <- lm_imp(SBP ~ age + gender + log(bili) + exp(creat),
                models = c(bili = 'glm_gamma_log', creat = 'lm'), data = NHANES)

If only one-sided truncation is needed, the other limit can be provided as NA.

Functions that are not available in R

It is possible to use functions that have different names in R and JAGS, or that do exist in JAGS, but not in R, by defining a new function in R that has the name of the function in JAGS.

Example:
In JAGS the inverse logit transformation is defined in the function ilogit. In R, there is no function ilogit, but the inverse logit is available as the distribution function of the logistic distribution plogis().

# Define the function ilogit
ilogit <- plogis

# Use ilogit in the model formula
mod5a <- lm_imp(SBP ~ age + gender + ilogit(creat), data = NHANES)

Nested functions

It is also possible to nest a function in another function.

Example:²

The complementary log log transformation is restricted to values larger than 0 and smaller than 1. In order to use this function on a variable that exceeds this range, as is the case for creat, a second transformation might be used, for instance the inverse logit from the previous example.

In JAGS, the complementary log-log transformation is implemented as cloglog, but since this function does not exist in (base) R, we first need to define it:

# define the complementary log log transformation
cloglog <- function(x) log(-log(1 - x))

# define the inverse logit (in case you have not done it in the previous example)
ilogit <- plogis

# nest ilogit inside cloglog
mod6a <- lm_imp(SBP ~ age + gender + cloglog(ilogit(creat)), data = NHANES)

Why do we need models for completely observed covariates?

The joint distribution of an outcome y, covariates x, random effects b and parameters θ, p(y, x, b, θ), is modelled as the product of univariate conditional distributions. To facilitate the specification of these distributions they are ordered so that longitudinal (level-1) variables may have baseline (level-2) variables in their linear predictors but not vice versa.

For example: $$\begin{align} p(y, x, b, \theta) = & p(y \mid x_1, ..., x_4, b_y, \theta_y) && \text{analysis model}\\ & p(x_1\mid \theta_{x1}) && \text{model for a complete baseline covariate}\\ & p(x_2\mid x_1, \theta_{x2}) && \text{model for an incomplete baseline covariate}\\ & p(x_3\mid x_1, x_2, b_{x3}, \theta_{x3}) && \text{model for a complete longitudinal covariate}\\ & p(x_4\mid x_1, x_2, x_3, b_{x4}, \theta_{x4}) && \text{model for an incomplete longitudinal covariate}\\ & p(b_y|\theta_b) p(b_{x3}|\theta_b) p(b_{x4}|\theta_b) && \text{models for the random effects}\\ & p(\theta_y) p(\theta_{x1}) \ldots p(\theta_{x4}) p(\theta_b) && \text{prior distributions}\end{align}$$

Since the parameter vectors θ_x1, θ_x2, … are assumed to be a priori independent, and furthermore x₁ is completely observed and modelled independently of incomplete variables, estimation of the other model parts is not affected by p(x₁ ∣ θ_x1) and, hence, this model can be omitted.

p(x₃ ∣ x₁, x₂, b_x3, θ_x3), on the other hand is modelled conditional on the incomplete covariate x₂ and can therefore not be omitted.

If there were no incomplete baseline covariates, i.e., if x₂ was completely observed, p(x₃ ∣ x₁, x₂, b_x3, θ_x3) would also not affect the estimation of parameters in the other parts of the model and could be omitted.

Setting reference categories for all variables

To specify the choice of reference category globally for all variables in the model, refcats can be set as

• refcats = "first"
• refcats = "last"
• refcats = "largest"

For example:

mod10a <- lm_imp(SBP ~ gender + age + race + educ + occup + smoke,
                 refcats = "largest", data = NHANES, n.adapt = 0)
#> Warning: 
#> It is currently not possible to use "contr.poly" for incomplete categorical covariates. I
#> will use "contr.treatment" instead.  You can specify (globally) which types of contrasts
#> are used by changing "options('contrasts')".

Setting reference categories for individual variables

Alternatively, refcats takes a named vector, in which the reference category for each variable can be specified either by its number or its name, or one of the three global types: “first”, “last” or “largest”. For variables for which no reference category is specified in the list the default is used.

mod10b <- lm_imp(SBP ~ gender + age + race + educ + occup + smoke,
                 refcats = list(occup = "not working", race = 3, educ = 'largest'),
                 data = NHANES, n.adapt = 0)
#> Warning: 
#> It is currently not possible to use "contr.poly" for incomplete categorical covariates. I
#> will use "contr.treatment" instead.  You can specify (globally) which types of contrasts
#> are used by changing "options('contrasts')".

To help to specify the reference category, the function set_refcat() can be used. It prints the names of the categorical variables that are selected by

• a specified model formula and/or
• a vector of auxiliary variables, or
• a vector of naming covariates

or all categorical variables in the data if only data is provided, and asks a number of questions which the user needs to reply to by input of a number.

refs_mod10 <- set_refcat(NHANES, formula = formula(mod10b))
#> The categorical variables are:
#> - "gender"
#> - "race"
#> - "educ"
#> - "occup"
#> - "smoke"
#> 
#> How do you want to specify the reference categories?
#> 
#> 1: Use the first category for each variable.
#> 2: Use the last category for each variabe.
#> 3: Use the largest category for each variable.
#> 4: Specify the reference categories individually.

When option 4 is chosen, a question for each categorical variable is asked, for example:

#> The reference category for “race” should be
#>
#> 1: Mexican American
#> 2: Other Hispanic
#> 3: Non-Hispanic White
#> 4: Non-Hispanic Black
#> 5: other

After specification of the reference categories for all categorical variables, the determined specification for the argument refcats is printed:

#> In the JointAI model specify:
#>  refcats = c(gender = 'female', race = 'Non-Hispanic White', educ = 'low',
#>              occup = 'not working', smoke = 'never')
#>
#> or use the output of this function.

set_refcat() also returns a named vector that can be passed to the argument refcats:

mod10c <- lm_imp(SBP ~ gender + age + race + educ + occup + smoke,
                 refcats = refs_mod10, data = NHANES, n.adapt = 0)
#> Warning: 
#> It is currently not possible to use "contr.poly" for incomplete categorical covariates. I
#> will use "contr.treatment" instead.  You can specify (globally) which types of contrasts
#> are used by changing "options('contrasts')".

Note:
Changing a reference category via the argument refcats does not change the order of levels in the dataset or any of the data matrices inside JointAI. Only when, in the JAGS model, the categorical variables is converted into dummy variables, the reference category is used to determine for which levels the dummies are created.