User-Defined Family Objects

In addition to the standard family objects in R, function mixed_model() also allows for user-specified family objects to be given in its family argument. The minimum requirements are that user-specified family object is of class "family" and is a list with the following components:

  • family: a character string giving the name of the family.

  • link: a character string giving the name of the link function.

  • linkfun: a function that transforms the mean of the distribution of the outcome to the linear predictor scale.

  • linkinv: a function that is the inverse of the linkfun defined above.

  • log_dens: a function that returns the log of the probability density function or the log of the probability mass function of the distribution of the outcome given the random effects. This function should have exactly five arguments named, y the repeated measurements outcome, eta denoting the linear predictor (fixed + random effects), mu_fun the function to calculate the mean from the linear predictor (i.e., the linkinv function), phis a vector of possible extra dispersion/scale parameters, and eta_zi the linear predictor (fixed + random effects) for the logistic regression for the extra zeros; relevant in the cases of zero-inflated and two-part models. Even if phis and eta_zi may not be required, it should always appear as arguments of log_dens. Moreover, when the extra parameters phis are required, they should be appropriately transformed in order for the optimization of the log-likelihood to be unconstrained with regard to these parameters. Finally, log_dens should be fully vectorized with respect to eta and eta_zi. That is, for a numeric vector y and a numeric matrix eta, it should return a matrix with the log density evaluated at each column of eta, and likewise for eta_zi.

For the internal calculations the derivative of log_dens with respect to the linear predictor eta and extra parameters phis (when present) is also required. If it is not given, then it is approximated numerically using a central difference approximation. However, it greatly facilitates stability and speed if the user also specifies these functions. If the user specifies these functions, they should have the names score_eta_fun and score_phis_fun, respectively, and both of them have three arguments named y the numeric vector of the repeated measurements outcome, mu the numeric vector or matrix of the mean of the outcome and phis a numeric vector of the extra parameters.

Negative Binomial Mixed Effects Model

We first illustrate the use of a custom-made family object in the case of a negative binomial repeated measurements outcome. We should note however that a family object negative.binomial() is already defined in the package and users could invoke this instead - for more info see the relative online help page for negative.binomial(). We start by simulating some data:

set.seed(101)
dd <- expand.grid(f1 = factor(1:3), f2 = LETTERS[1:2], g = 1:30, rep = 1:15,
                  KEEP.OUT.ATTRS = FALSE)
mu <- 5*(-4 + with(dd, as.integer(f1) + 4 * as.numeric(f2)))
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)

Next we define a new family object for negative binomial data

And we fit the model

In the call to mixed_model() we specify the argument n_phis which lets the function know how many extra parameters does the distribution of the outcome have. In addition, we also specify that better quality initial values for the fixed effects betas can be obtained by fitting a Poisson regression model to the repeated measurements outcome (ignoring the correlations) with the same design matrix as the fixed-effects design matrix. For the unique elements of the random effects covariance matrix and for the extra parameters phis the default initial values are used.

In the output, the reported phi parameters is on the log scale, because in the definition of log_dens we specified that phis <- exp(phis).

As mentioned above, the internal computations are facilitated if the user also specifies the derivatives of the log density with respect to the linear predictor eta and the extra parameters phis. Hence, in the following we define a new negative binomial family object with these functions also defined:

And we fit the model

We obtain almost identical results, but in half of the computing time in this example. Depending on the size of the data set, the gain in computing time when the derivatives for eta and phis are provided can be much more substantial.

Student’s t Mixed Effects Model

We further illustrate the use of custom-made family objects by fitting a linear mixed model but with Student’s t distributed error terms. We start by simulating some data

Following we define the custom-made family object for the Student’s t distribution:

And we fit the model (results not shown)

fm <-  mixed_model(y ~ sex * time, random = ~ time | id, data = DF, 
                   family = students.t(4), n_phis = 1, 
                   initial_values = list("betas" = gaussian()))

As was the case also for the negative binomial model, the extra phis parameter has been appropriately transformed to have an unconstrained optimization of the likelihood. Hence, the estimated sigma parameter is given by exp(fm$phis).

Beta Mixed-Effects Model

As a final illustration we show how a beta mixed effects model can be fitted, i.e., applicable in case of bounded repeated measurements outcomes. We start again by simulating some data:

Next we define the custom-made family object as in the case of the Student’s t distribution:

And we fit the model

To depict the results of the model we create an effects plot, showing the longitudinal evolution of the average/median male and female. We start by constructing the data frame that contains the values we want to depict, and using it in the effectPlotData() function; just for illustration, sandwich/robust standard errors are used in the calculation of the 95% poinwise cofidence intervals:

nDF <- with(DF, expand.grid(time = seq(min(time), max(time), length = 25),
                            sex = levels(sex)))

plot_data <- effectPlotData(gm, nDF, sandwich = TRUE)

The figure is created with a call to xyplot() from the lattice package: