glmmrBase (version 0.2.5)

Model: A GLMM Model

Description

A GLMM Model

A GLMM Model

Public fields

covariance

A Covariance object defining the random effects covariance.

mean_function

A MeanFunction object, defining the mean function for the model, including the data and covariate design matrix X.

family

One of the family function used in R's glm functions. See family for details

exp_condition

A vector indicting the unique experimental conditions for each observation, see Details.

Sigma

The overall covariance matrix for the observations. Calculated and updated automatically as \(W^{-1} + ZDZ^T\) where W is an n x n diagonal matrix with elements on the diagonal equal to the GLM iterated weights. See Details.

var_par

Scale parameter required for some distributions (Gaussian, Gamma, Beta).

Methods

Public methods

• Model$use_attenuation()

• Model$fitted()

• Model$new()

• Model$print()

• Model$n()

• Model$n_cluster()

• Model$subset_rows()

• Model$subset_cols()

• Model$sim_data()

• Model$check()

• Model$update_parameters()

• Model$information_matrix()

• Model$power()

• Model$clone()

---

Method use_attenuation()

Sets the model to use or not use "attenuation" when calculating the first-order approximation to the covariance matrix.

Usage

Model$use_attenuation(use)

Arguments

use

Logical indicating whether to use "attenuation".

Returns

None. Used for effects.

---

Method fitted()

Return predicted values based on the currently stored parameter values in `mean_function`

Usage

Model$fitted(type = "link")

Arguments

type

One of either "`link`" for values on the scale of the link function, or "`response`" for values on the scale of the response

Returns

A Matrix class object containing the predicted values

---

Method new()

Create a new Model object

Usage

Model$new(
  covariance,
  mean,
  data = NULL,
  family = NULL,
  var_par = NULL,
  verbose = TRUE,
  skip.sigma = FALSE
)

Arguments

covariance

Either a Covariance object, or an equivalent list of arguments that can be passed to `Covariance` to create a new object. At a minimum the list must specify a formula. If parameters are not included then they are initialised to 0.5.

mean

Either a MeanFunction object, or an equivalent list of arguments that can be passed to `MeanFunction` to create a new object. At a minimum the list must specify a formula. If parameters are not included then they are initialised to 0.

data

A data frame with the data required for the mean function and covariance objects. This argument can be ignored if data are provided to the covariance or mean arguments either via `Covariance` and `MeanFunction` object, or as a member of the list of arguments to both `covariance` and `mean`.

family

A family object expressing the distribution and link function of the model, see family. This argument is optional if the family is provided either via a `MeanFunction` or `MeanFunction` objects, or as members of the list of arguments to `mean`. Current accepts binomial, gaussian, Gamma, poisson, and Beta.

var_par

Scale parameter required for some distributions, including Gaussian. Default is NULL.

verbose

Logical indicating whether to provide detailed output

skip.sigma

Logical indicating whether to skip the creating of the covariance matrix Sigma. For very large designs with thousands of observations or more, the covariance matrix will be too big to fit in memory, so this option will prevent sigma being created.

Returns

A new Model class object

Examples

#create a data frame describing a cross-sectional parallel cluster
#randomised trial
df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
mf1 <- MeanFunction$new(
  formula = ~ factor(t) + int - 1,
  data=df
)
cov1 <- Covariance$new(
  data = df,
  formula = ~ (1|gr(cl)) + (1|gr(cl*t))
)
des <- Model$new(
  covariance = cov1,
  mean = mf1,
  family = stats::gaussian(),
  var_par = 1
)
#alternatively we can pass the data directly to Model
#here we will specify a cohort study and provide parameter values
df <- nelder(~ind(20) * t(6))
df$int <- 0
df[df$t > 3, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(ind)),
    parameters = c(0.25)),
  mean = list(
    formula = ~ int,
    parameters = c(1,0.5)),
  data = df,
  family = stats::poisson())
                  
#an example of a spatial grid with two time points
df <- nelder(~ (x(10)*y(10))*t(2))
spt_design <- Model$new(covariance = list( formula = ~(1|fexp(x,y)*ar1(t))),
                         mean = list(formula = ~ 1),
                         data = df,
                         family = stats::gaussian()) 

---

Method print()

Print method for `Model` class

Usage

Model$print()

Arguments

...

ignored

---

Method n()

Returns the number of observations in the model

Usage

Model$n(...)

Arguments

...

ignored

---

Method n_cluster()

Returns the number of clusters at each level

Usage

Model$n_cluster()

Arguments

...

ignored

Returns

A data frame with the level, number of clusters, and variables describing each level.

Examples

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(formula = ~ (1|gr(cl)) + (1|gr(cl*t))),
  mean = list(formula = ~ factor(t) + int - 1),
  data = df, 
  family = stats::gaussian(),
  var_par = 1
)
des$n_cluster() ## returns two levels of 10 and 50

---

Method subset_rows()

Subsets the design keeping specified observations only

Given a vector of row indices, the corresponding rows will be kept and the other rows will be removed from the mean function and covariance

Usage

Model$subset_rows(index)

Arguments

index

Integer or vector integers listing the rows to keep

Returns

The function updates the object and nothing is returned

Examples

#generate a stepped wedge design and remove the first sequence
des <- stepped_wedge(8,10,icc=0.05)
ids_to_keep <- which(des$mean_function$data$J!=1)
des$subset_rows(ids_to_keep)

---

Method subset_cols()

Subsets the columns of the design

Removes the specified columns from the linked mean function object's X matrix.

Usage

Model$subset_cols(index)

Arguments

index

Integer or vector of integers specifying the indexes of the columns to keep

Returns

The function updates the object and nothing is returned

Examples

#generate a stepped wedge design and remove first and last time periods
des <- stepped_wedge(8,10,icc=0.05)
des$subset_cols(c(2:8,10))

---

Method sim_data()

Generates a realisation of the design

Generates a single vector of outcome data based upon the specified GLMM design

Usage

Model$sim_data(type = "y")

Arguments

type

Either 'y' to return just the outcome data, 'data' to return a data frame with the simulated outcome data alongside the model data, or 'all', which will return a list with simulated outcomes y, matrices X and Z, parameters beta, and the values of the simulated random effects.

Returns

Either a vector, a data frame, or a list

Examples

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(cl)*ar1(t)),
    parameters = c(0.25,0.8)),
  mean = list(
    formula = ~ factor(t) + int - 1,
    parameters = c(rep(0,5),0.6)),
  data = df,
  family = stats::binomial()
)
ysim <- des$sim_data()

---

Method check()

Checks for any changes in linked objects and updates.

Checks for any changes in any object and updates all linked objects if any are detected. Generally called automatically.

Usage

Model$check(verbose = TRUE)

Arguments

verbose

Logical indicating whether to report if any updates are made, defaults to TRUE

Returns

Linked objects are updated by nothing is returned

Examples

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(cl)*ar1(t)),
    parameters = c(0.25,0.8)),
  mean = list(
    formula = ~ factor(t) + int - 1,
    parameters = c(rep(0,5),0.6)),
  data = df, 
  family = stats::binomial()
)
des$check() #does nothing
des$covariance$parameters <- c(0.1,0.9)
des$check() #updates 
des$mean_function$parameters <- c(rnorm(5),0.1)
des$check() #updates 

---

Method update_parameters()

Updates the parameters of the mean function and/or the covariance function

Usage

Model$update_parameters(mean.pars = NULL, cov.pars = NULL, verbose = FALSE)

Arguments

mean.pars

(Optional) Vector of new mean function parameters

cov.pars

(Optional) Vector of new covariance function(s) parameters

verbose

Logical indicating whether to provide more detailed feedback

Examples

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(cl)*ar1(t))),
  mean = list(
    formula = ~ factor(t) + int - 1),
  data = df, 
  family = stats::binomial()
)
des$update_parameters(cov.pars = c(0.1,0.9))

---

Method information_matrix()

Generates the information matrix

Usage

Model$information_matrix()

Returns

A PxP matrix

---

Method power()

Estimates the power of the design described by the model using the square root of the relevant element of the GLS variance matrix:

$$(X^T\Sigma^{-1}X)^{-1}$$

Note that this is equivalent to using the "design effect" for many models.

Usage

Model$power(alpha = 0.05, two.sided = TRUE, alternative = "pos")

Arguments

alpha

Numeric between zero and one indicating the type I error rate. Default of 0.05.

two.sided

Logical indicating whether to use a two sided test

alternative

For a one-sided test whether the alternative hypothesis is that the parameter is positive "pos" or negative "neg"

Returns

A data frame describing the parameters, their values, expected standard errors and estimated power.

Examples

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(cl)) + (1|gr(cl*t)),
    parameters = c(0.25,0.1)),
  mean = list(
    formula = ~ factor(t) + int - 1,
    parameters = c(rep(0,5),0.6)),
  data = df, 
  family = stats::gaussian(),
  var_par = 1
)
des$power() #power of 0.90 for the int parameter

---

Method clone()

The objects of this class are cloneable with this method.

Usage

Model$clone(deep = FALSE)

Arguments

deep

Whether to make a deep clone.

Details

An R6 class representing a GLMM model

For the generalised linear mixed model

$$Y \sim F(\mu,\sigma)$$ $$\mu = h^-1(X\beta + Zu)$$ $$u \sim MVN(0,D)$$

where h is the link function. A Model in comprised of a MeanFunction object, which defines the family F, link function h, and fixed effects design matrix X, and a Covariance object, which defines Z and D. The class provides methods for analysis and simulation with these models.

This class provides methods for generating the matrices described above and data simulation, and serves as a base for extended functionality in related packages.

Many calculations use the covariance matrix of the observations, such as the information matrix, which is used in power calculations and other functions. For non-Gaussian models, the class uses the approximation proposed by Breslow and Clayton (1993) based on the marginal quasilikelihood:

$$\Sigma = W^{-1} + ZDZ^T$$

where _W_ is a diagonal matrix with the GLM iterated weights for each observation equal to, for individual _i_ \(\left( \frac{(\partial h^{-1}(\eta_i))}{\partial \eta_i}\right) ^2 Var(y|u)\) (see Table 2.1 in McCullagh and Nelder (1989) <ISBN:9780412317606>). For very large designs, this can be disabled as the memory requirements can be prohibitive (use option `skip_sigma`).

See glmmrBase for a detailed guide on model specification.

**Attenuation** For calculations such as the information matrix, the first-order approximation to the covariance matrix proposed by Breslow and Clayton (1993), described above, is used. The approximation is based on the marginal quasilikelihood. Zegers, Liang, and Albert (1988) suggest that a better approximation to the marginal mean is achieved by "attenuating" the linear predictor. Setting `use` equal to TRUE uses this adjustment for calculations using the covariance matrix for non-linear models.

Calls the respective print methods of the linked covariance and mean function objects.

The matrices X and Z both have n rows, where n is the number of observations in the model/design.

**Number of clusters** Returns a data frame describing the number of independent clusters or groups at each level in the design. For example, if there were cluster-periods nested in clusters, then the top level would be clusters, and the second level would be cluster periods.

Using `update_parameters()` is the preferred way of updating the parameters of the mean or covariance objects as opposed to direct assignment, e.g. `self$covariance$parameters <- c(...)`. The function calls check functions to automatically update linked matrices with the new parameters. If using direct assignment, call `self$check()` afterwards.

References

Breslow, N. E., Clayton, D. G. (1993). Approximate Inference in Generalized Linear Mixed Models. Journal of the American Statistical Association<, 88(421), 9–25. <doi:10.1080/01621459.1993.10594284>

Zeger, S. L., Liang, K.-Y., Albert, P. S. (1988). Models for Longitudinal Data: A Generalized Estimating Equation Approach. Biometrics, 44(4), 1049.<doi:10.2307/2531734>

See Also

nelder, MeanFunction, Covariance

Examples

## ------------------------------------------------
## Method `Model$new`
## ------------------------------------------------

#create a data frame describing a cross-sectional parallel cluster
#randomised trial
df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1

mf1 <- MeanFunction$new(
  formula = ~ factor(t) + int - 1,
  data=df
)
cov1 <- Covariance$new(
  data = df,
  formula = ~ (1|gr(cl)) + (1|gr(cl*t))
)
des <- Model$new(
  covariance = cov1,
  mean = mf1,
  family = stats::gaussian(),
  var_par = 1
)

#alternatively we can pass the data directly to Model
#here we will specify a cohort study and provide parameter values
df <- nelder(~ind(20) * t(6))
df$int <- 0
df[df$t > 3, 'int'] <- 1

des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(ind)),
    parameters = c(0.25)),
  mean = list(
    formula = ~ int,
    parameters = c(1,0.5)),
  data = df,
  family = stats::poisson())
                  
#an example of a spatial grid with two time points
df <- nelder(~ (x(10)*y(10))*t(2))
spt_design <- Model$new(covariance = list( formula = ~(1|fexp(x,y)*ar1(t))),
                         mean = list(formula = ~ 1),
                         data = df,
                         family = stats::gaussian()) 

## ------------------------------------------------
## Method `Model$n_cluster`
## ------------------------------------------------

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1

des <- Model$new(
  covariance = list(formula = ~ (1|gr(cl)) + (1|gr(cl*t))),
  mean = list(formula = ~ factor(t) + int - 1),
  data = df, 
  family = stats::gaussian(),
  var_par = 1
)
des$n_cluster() ## returns two levels of 10 and 50

## ------------------------------------------------
## Method `Model$subset_rows`
## ------------------------------------------------

#generate a stepped wedge design and remove the first sequence
des <- stepped_wedge(8,10,icc=0.05)
ids_to_keep <- which(des$mean_function$data$J!=1)
des$subset_rows(ids_to_keep)

## ------------------------------------------------
## Method `Model$subset_cols`
## ------------------------------------------------

#generate a stepped wedge design and remove first and last time periods
des <- stepped_wedge(8,10,icc=0.05)
des$subset_cols(c(2:8,10))

## ------------------------------------------------
## Method `Model$sim_data`
## ------------------------------------------------

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(cl)*ar1(t)),
    parameters = c(0.25,0.8)),
  mean = list(
    formula = ~ factor(t) + int - 1,
    parameters = c(rep(0,5),0.6)),
  data = df,
  family = stats::binomial()
)
ysim <- des$sim_data()

## ------------------------------------------------
## Method `Model$check`
## ------------------------------------------------

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(cl)*ar1(t)),
    parameters = c(0.25,0.8)),
  mean = list(
    formula = ~ factor(t) + int - 1,
    parameters = c(rep(0,5),0.6)),
  data = df, 
  family = stats::binomial()
)
des$check() #does nothing
des$covariance$parameters <- c(0.1,0.9)
des$check() #updates 
des$mean_function$parameters <- c(rnorm(5),0.1)
des$check() #updates 

## ------------------------------------------------
## Method `Model$update_parameters`
## ------------------------------------------------

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(cl)*ar1(t))),
  mean = list(
    formula = ~ factor(t) + int - 1),
  data = df, 
  family = stats::binomial()
)
des$update_parameters(cov.pars = c(0.1,0.9))

## ------------------------------------------------
## Method `Model$power`
## ------------------------------------------------

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    formula = ~ (1|gr(cl)) + (1|gr(cl*t)),
    parameters = c(0.25,0.1)),
  mean = list(
    formula = ~ factor(t) + int - 1,
    parameters = c(rep(0,5),0.6)),
  data = df, 
  family = stats::gaussian(),
  var_par = 1
)
des$power() #power of 0.90 for the int parameter