Learn R Programming

glmmrBase (version 0.1.2)

Model: A GLMM Model

Description

A GLMM Model

A GLMM Model

Arguments

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.

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


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.function,
  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.

mean.function

Either a MeanFunction object, or an equivalent list of arguments that can be passed to `MeanFunction` to create a new object.

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, parameters = c(rep(0,5),0.6), family = stats::gaussian() ) cov1 <- Covariance$new( data = df, formula = ~ (1|gr(cl)) + (1|gr(cl*t)), parameters = c(0.25,0.1) ) des <- Model$new( covariance = cov1, mean.function = mf1, var_par = 1 )

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

des <- Model$new( covariance = list( data=df, formula = ~ (1|gr(ind)*ar1(t)), parameters = c(0.25,0.8)), mean.function = list( formula = ~factor(t) + int - 1, data=df, parameters = rep(0,7), 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(data=df, formula = ~(1|fexp(x,y)*ar1(t)), parameters =c(0.2,0.1,0.8)), mean.function = list(data=df, formula = ~ 1, parameters = c(0.5), family=stats::poisson()))


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

mf1 <- MeanFunction$new( formula = ~ factor(t) + int - 1, data=df, parameters = c(rep(0,5),0.6), family = stats::gaussian() ) cov1 <- Covariance$new( data = df, formula = ~ (1|gr(cl)) + (1|gr(cl*t)), parameters = c(0.25,0.1) ) des <- Model$new( covariance = cov1, mean.function = mf1, 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, or 'data' to return a data frame with the simulated outcome data alongside the model data

Returns

Either a vector or a data frame

Examples

df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
  covariance = list(
    data = df,
    formula = ~ (1|gr(cl)*ar1(t)),
    parameters = c(0.25,0.8)),
  mean.function = list(
    formula = ~ factor(t) + int - 1,
    data=df,
    parameters = c(rep(0,5),0.6),
    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(
    data = df,
    formula = ~ (1|gr(cl)*ar1(t)),
    parameters = c(0.25,0.8)),
  mean.function = list(
    formula = ~ factor(t) + int - 1,
    data=df,
    parameters = c(rep(0,5),0.6),
    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 information_matrix()

Generates the information matrix

Usage

Model$information_matrix()

Returns

A PxP matrix


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 + Z\gamma)$$ $$\gamma \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.

The class by default calculates the covariance matrix of the observations as:

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

where _W_ is a diagonal matrix with the WLS iterated weights for each observation equal to, for individual _i_ \(\phi a_i v(\mu_i)[h'(\mu_i)]^2\) (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.

See glmmrBase for a detailed guide on model specification.

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.

See Also

nelder, MeanFunction, Covariance

Examples

Run this code

## ------------------------------------------------
## 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,
  parameters = c(rep(0,5),0.6),
  family = stats::gaussian()
)
cov1 <- Covariance$new(
  data = df,
  formula = ~ (1|gr(cl)) + (1|gr(cl*t)),
  parameters = c(0.25,0.1)
)
des <- Model$new(
  covariance = cov1,
  mean.function = mf1,
  var_par = 1
)

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

des <- Model$new(
covariance = list(
  data=df,
  formula = ~ (1|gr(ind)*ar1(t)),
  parameters = c(0.25,0.8)),
mean.function = list(
  formula = ~factor(t) + int - 1,
  data=df,
  parameters = rep(0,7),
  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(data=df,
                                           formula = ~(1|fexp(x,y)*ar1(t)),
                                           parameters =c(0.2,0.1,0.8)),
                         mean.function = list(data=df,
                                              formula = ~ 1,
                                              parameters = c(0.5),
                                              family=stats::poisson())) 

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

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,
  parameters = c(rep(0,5),0.6),
  family = stats::gaussian()
)
cov1 <- Covariance$new(
  data = df,
  formula = ~ (1|gr(cl)) + (1|gr(cl*t)),
  parameters = c(0.25,0.1)
)
des <- Model$new(
  covariance = cov1,
  mean.function = mf1,
  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(
    data = df,
    formula = ~ (1|gr(cl)*ar1(t)),
    parameters = c(0.25,0.8)),
  mean.function = list(
    formula = ~ factor(t) + int - 1,
    data=df,
    parameters = c(rep(0,5),0.6),
    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(
    data = df,
    formula = ~ (1|gr(cl)*ar1(t)),
    parameters = c(0.25,0.8)),
  mean.function = list(
    formula = ~ factor(t) + int - 1,
    data=df,
    parameters = c(rep(0,5),0.6),
    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 

Run the code above in your browser using DataLab