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
useLogical indicating whether to use "attenuation".
Returns
None. Used for effects.
Return fitted values. Does not account for the random effects. For simulated values based
on resampling random effects, see also sim_data(). To predict the values including random effects at a new location see also
predict().
Usage
Model$fitted(type = "link", X, u, sample = FALSE, sample_n = 100)
Arguments
typeOne of either "link" for values on the scale of the link function, or "response"
for values on the scale of the response
X
(Optional) Fixed effects matrix to generate fitted values
u
(Optional) Random effects values at which to generate fitted values
sample
Logical. If TRUE then the parameters will be re-sampled from their sampling distribution. Currently only works
with existing X matrix and not user supplied matrix X and this will also ignore any provided random effects.
sample_n
Integer. If sample is TRUE, then this is the number of samples.
Returns
Fitted values as either a vector or matrix depending on the number of samples
Generates the residuals for the model
Generates one of several types of residual for the model. If conditional = TRUE then
the residuals include the random effects, otherwise only the fixed effects are included. For type,
there are raw, pearson, and standardized residuals. For conditional residuals a matrix is returned
with each column corresponding to a sample of the random effects.
Usage
Model$residuals(type = "standardized", conditional = TRUE)
Arguments
typeEither "standardized", "raw" or "pearson"
conditional
Logical indicating whether to condition on the random effects (TRUE) or not (FALSE)
Returns
A matrix with either one column is conditional is false, or with number of columns corresponding
to the number of MCMC samples.
Generate predictions at new values
Generates predicted values using a new data set to specify covariance
values and values for the variables that define the covariance function.
The function will return a list with the linear predictor, conditional
distribution of the new random effects term conditional on the current estimates
of the random effects, and some simulated values of the random effects if requested.
Usage
Model$predict(newdata, offset = rep(0, nrow(newdata)), m = 0)
Arguments
newdataA data frame specifying the new data at which to generate predictions
offset
Optional vector of offset values for the new data
m
Number of samples of the random effects to draw
Returns
A list with the linear predictor, parameters (mean and covariance matrices) for
the conditional distribution of the random effects, and any random effect samples.
Method new()
Create a new Model object. Typically, a model is generated from a formula and data. However, it can also be
generated from a previous model fit.
Usage
Model$new(
formula,
covariance,
mean,
data = NULL,
family = NULL,
var_par = NULL,
offset = NULL,
weights = NULL,
trials = NULL,
model_fit = NULL
)
Arguments
formulaA model formula containing fixed and random effect terms. The formula can be one way (e.g. ~ x + (1|gr(cl))) or
two-way (e.g. y ~ x + (1|gr(cl))). One way formulae will generate a valid model enabling data simulation, matrix calculation,
and power, etc. Outcome data can be passed directly to model fitting functions, or updated later using member function update_y().
For binomial models, either the syntax cbind(y, n-y) can be used for outcomes, or just y and the number of trials passed to the argument
trials described below.
covariance
(Optional) Either a Covariance object, an equivalent list of arguments
that can be passed to Covariance to create a new object, or a vector of parameter values. At a minimum the list must specify a formula.
If parameters are not included then they are initialised to 0.5.
mean
(Optional) Either a MeanFunction object, an equivalent list of arguments
that can be passed to MeanFunction to create a new object, or a vector of parameter values. 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. Currently accepts binomial,
gaussian, Gamma, poisson, Beta, and Quantile.
var_par
(Optional) Scale parameter required for some distributions, including Gaussian. Default is NULL.
offset
(Optional) A vector of offset values. Optional - could be provided to the argument to mean instead.
weights
(Optional) A vector of weights.
trials
(Optional) For binomial family models, the number of trials for each observation. If it is not set, then it will
default to 1 (a bernoulli model).
model_fit
(optional) A mcml model fit resulting from a call to MCML or LA
Returns
A new Model class object
Examples
\dontshow{
setParallel(FALSE)
}
# For more examples, see the examples for MCML.#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
mod <- Model$new(
formula = ~ factor(t) + int - 1 + (1|gr(cl)) + (1|gr(cl,t)),
data = df,
family = stats::gaussian()
)
# We can also include the outcome data in the model initialisation.
# For example, simulating data and creating a new object:
df$y <- mod$sim_data()
mod <- Model$new(
formula = y ~ factor(t) + int - 1 + (1|gr(cl)) + (1|gr(cl,t)),
data = df,
family = stats::gaussian()
)
# 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(
formula = ~ int + (1|gr(ind)),
data = df,
family = stats::poisson()
)
# or with parameter values specified
des <- Model$new(
formula = ~ int + (1|gr(ind)),
covariance = c(0.05),
mean = 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(formula = ~ 1 + (1|ar0(t)*fexp(x,y)),
data = df,
family = stats::gaussian())
Print method for Model class
Method n()
Returns the number of observations in the model
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
indexInteger or vector integers listing the rows to keep
Returns
The function updates the object and nothing is returned.
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
typeEither '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
\dontshow{
setParallel(FALSE) # for the CRAN check
}
des <- Model$new(
formula = ~ factor(t) + int - 1 + (1|gr(cl)*ar0(t)),
covariance = c(0.05,0.8),
mean = c(rep(0,5),0.6),
data = df,
family = stats::binomial()
)
ysim <- des$sim_data()
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, var.par = NULL)
Arguments
mean.pars(Optional) Vector of new mean function parameters
cov.pars
(Optional) Vector of new covariance function(s) parameters
var.par
(Optional) A scalar value for var_par
Examples
\dontshow{
setParallel(FALSE) # for the CRAN check
}
df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
formula = ~ factor(t) + int - 1 + (1|gr(cl)*ar0(t)),
data = df,
family = stats::binomial()
)
des$update_parameters(cov.pars = c(0.1,0.9))
Arguments
include.relogical indicating whether to return the information matrix including the random effects components (TRUE),
or the mixed model information matrix for beta only (FALSE).
theta
Logical. If TRUE the function will return the variance-coviariance matrix for the covariance parameters and ignore the first argument. Otherwise, the fixed effect
parameter information matrix is returned.
oim
Logical. If TRUE, returns the observed information matrix for both beta and theta, disregarding other arguments to the function.
Method sandwich()
Returns the robust sandwich variance-covariance matrix for the fixed effect parameters
Method small_sample_correction()
Returns a small sample correction. The option "KR" returns the Kenward-Roger bias-corrected variance-covariance matrix
for the fixed effect parameters and degrees of freedom. Option "KR2" returns an improved correction given
in Kenward & Roger (2009) doi:j.csda.2008.12.013. Note, that the corrected/improved version is invariant
under reparameterisation of the covariance, and it will also make no difference if the covariance is linear
in parameters. Exchangeable covariance structures in this package (i.e. gr()) are parameterised in terms of
the variance rather than standard deviation, so the results will be unaffected. Option "sat" returns the "Satterthwaite"
correction, which only includes corrected degrees of freedom, along with the GLS standard errors.
Usage
Model$small_sample_correction(type, oim = FALSE)
Arguments
typeEither "KR", "KR2", or "sat", see description.
oim
Logical. If TRUE use the observed information matrix, otherwise use the expected information matrix
Returns the inferential statistics (F-stat, p-value) for a modified Box correction doi:10.1002/sim.4072 for
Gaussian-identity models.
Arguments
yOptional. If provided, will update the vector of outcome data. Otherwise it will use the data from
the previous model fit.
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
alphaNumeric 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
\dontshow{
setParallel(FALSE) # for the CRAN check
}
df <- nelder(~(cl(10)*t(5)) > ind(10))
df$int <- 0
df[df$cl > 5, 'int'] <- 1
des <- Model$new(
formula = ~ factor(t) + int - 1 + (1|gr(cl)) + (1|gr(cl,t)),
covariance = c(0.05,0.1),
mean = 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 w_matrix()
Returns the diagonal of the matrix W used to calculate the covariance matrix approximation
Returns
A vector with values of the glm iterated weights
Method dh_deta()
Returns the derivative of the link function with respect to the linear preditor
Method Sigma()
Returns the (approximate) covariance matrix of y
Returns the covariance matrix Sigma. For non-linear models this is an approximation. See Details.
Usage
Model$Sigma(inverse = FALSE)
Arguments
inverseLogical indicating whether to provide the covariance matrix or its inverse
Method MCML()
Stochastic Maximum Likelihood model fitting
Usage
Model$MCML(
y = NULL,
method = "saem",
tol = 0.01,
max.iter = 50,
se = "gls",
oim = FALSE,
reml = FALSE,
mcmc.pkg = "rstan",
se.theta = TRUE,
iter.warmup = NULL,
iter.sampling = NULL,
chains = NULL,
lower.bound = NULL,
upper.bound = NULL,
lower.bound.theta = NULL,
upper.bound.theta = NULL,
alpha = 0.8,
convergence.prob = 0.95,
pr.average = FALSE,
bf.tol = 10,
bf.hist = 10,
bf.k0 = 10,
importance = FALSE,
conv.criterion = 2,
skip.theta = FALSE,
constr.zero = 1
)
Arguments
yOptional. A numeric vector of outcome data. If this is not provided then either the outcome must have been specified when
initialising the Model object, or the outcome data has been updated using member function update_y()
method
The MCML algorithm to use, either mcem or mcnr, mcnr2, or saem see Details. Default is saem. mcem.adapt, mcnr2.adapt and mcnr.adapt will use adaptive
MCMC sample sizes starting small and increasing to the the maximum value specified in mcmc_options$sampling, which may result in faster convergence. saem uses a
stochastic approximation expectation maximisation algorithm. MCMC samples are kept from all iterations and so a smaller number of samples are needed per iteration. The
qualifier .dual can also be added (e.g. saem.dual), which combines the fixed and covariance parameter estimation steps.
tol
Numeric value, tolerance of the MCML algorithm, the maximum difference in parameter estimates
between iterations at which to stop the algorithm. If two values are provided then different tolerances will be
applied to the fixed effect and covariance parameters.
max.iter
Integer. The maximum number of iterations of the MCML algorithm.
se
String. Type of standard error and/or inferential statistics to return. Options are "gls" for GLS standard errors (the default),
"robust" for robust standard errors, "kr" for original Kenward-Roger bias corrected standard errors,
"kr2" for the improved Kenward-Roger correction, "sat" for Satterthwaite degrees of freedom correction (this is the same
degrees of freedom correction as Kenward-Roger, but with GLS standard errors), "box" to use a modified Box correction (does not return confidence intervals),
"bw" to use GLS standard errors with a between-within correction to the degrees of freedom, "bwrobust" to use robust
standard errors with between-within correction to the degrees of freedom.
oim
Logical. If TRUE use the observed information matrix, otherwise use the expected information matrix for standard error and related calculations.
reml
Logical. Whether to use a restricted maximum likelihood correction for fitting the covariance parameters
mcmc.pkg
String. Either rstan to use rstan sampler, or
analytic to use a Normal approximation to the posterior with direct estimates of the posterior mean and variance. cmdstanr will compile the MCMC programs to the library folder the first time they are run,
so may not currently be an option for some users.
se.theta
Logical. Whether to calculate the standard errors for the covariance parameters. This step is a slow part
of the calculation, so can be disabled if required in larger models. Has no effect for Kenward-Roger standard errors.
iter.warmup
Integer. The number of warmup iterations for each MCMC run on each iteration of the algorithm. If this value is left as NULL then the value stored in self$mcmc_options$warmup will be used.
iter.sampling
Integer. The number of sampling iterations for each MCMC run on each iteration of the algorithm. The default values have been selected to provide
relatively good convergence for the default SAEM algorithm, but may need to be increased for MCEM and MCNR. If an adaptive algorithm is used, then this is the maximum
number of iterations per MCMC run. If this value is left as NULL then the value stored in self$mcmc_options$samps will be used.
chains
Integer. The number of MCMC chains to run in parallel. The default is one, which generally provides good results. If this value is left as NULL then the value stored in self$mcmc_options$chains will be used.
lower.bound
Optional. Vector of lower bounds for the fixed effect parameters. To apply bounds use MCEM.
upper.bound
Optional. Vector of upper bounds for the fixed effect parameters. To apply bounds use MCEM.
lower.bound.theta
Optional. Vector of lower bounds for the covariance parameters (default is 0; negative values will cause an error)
upper.bound.theta
Optional. Vector of upper bounds for the covariance parameters.
alpha
If using SAEM then this parameter controls the step size. On each iteration i the step size is (1/alpha)^i, default is 0.8. Values around 0.5
will result in lower bias but slower convergence, values closer to 1 will result in higher convergence but potentially higher error.
convergence.prob
Numeric value in (0,1) indicating the probability of convergence if using convergence criteria 2, 3, or 4.
pr.average
Logical indicating whether to use Polyak-Ruppert averaging if using the SAEM algorithm (default is FALSE)
bf.tol
Integer indicating the Bayes Factor convergence criterion
bf.hist
Integer, the number of iterations in the running mean for the Bayes Factor convergence criterion
bf.k0
Integer, the expected number of iterations to convergence of the Bayes Factor convergence criterion.
importance
Logical. If TRUE and using the analytic samples, the u samples are weighted using an importance sampling step. If FALSE, it is equivalent to the Laplace
Approximation Gaussian distribution of the random effects.
conv.criterion
Integer. The convergence criterion for the algorithm. 1 = the maximum difference between parameter estimates between iterations as defined by tol,
2 = The probability of improvement in the overall log-likelihood is less than 1 - convergence.prob
3 = The probability of improvement in the log-likelihood for the fixed effects is less than 1 - convergence.prob
4 = The probabilities of improvement in the log-likelihood the fixed effects and covariance parameters are both less than 1 - convergence.prob
skip.theta
Logical. If TRUE then the covariance parameter estimation step is skipped. This option is mainly used for testing, but may be useful
if covariance parameters are known.
constr.zero
Scalar. A Soft sum-to-zero constraint can be forced on the random effects so that their sum is N(0,constr.zero*Q). Small values, e.g. 0.001
may be useful if there is possible identifiability issues for intercept terms, such as in more complex, or higher dimensional, random effects structures like spatial models.
Examples
\dontrun{
# Simulated trial data example
data(SimTrial,package = "glmmrBase")
model <- Model$new(
formula = y ~ int + factor(t) - 1 + (1|gr(cl)*ar1(t)),
data = SimTrial,
family = gaussian()
)
glm3 <- model$MCML()# Salamanders data example
data(Salamanders,package="glmmrBase")
model <- Model$new(
mating~fpop:mpop-1+(1|gr(mnum))+(1|gr(fnum)),
data = Salamanders,
family = binomial()
)
# use MCEM + REML with 500 sampling iterations
glm2 <- model$MCML(method = "mcem", iter.sampling = 500, reml = TRUE)
# as an alternative, we will specify the variance parameters on the
# log scale and use a fast fitting algorithm
# we will use two newton-raphson steps, and Normal approximation posteriors with
# conjugate gradient descent
# the maximum number of iterations is increased as it takes 100-110 in this example
# we can also chain together the functions
# specify starting values for the covariance parameters to prevent potential fit failure
# with random initialisation on the log scale
glm3 <- Model$new(
mating~fpop:mpop-1+(1|grlog(mnum))+(1|grlog(fnum)),
data = Salamanders,
family = binomial(),
covariance = c(-0.5,-0.5)
)$MCML(method = "mcnr2", mcmc.pkg = "analytic", iter.sampling = 50, max.iter = 150)
# Example using simulated data
#create example data with six clusters, five time periods, and five people per cluster-period
df <- nelder(~(cl(6)*t(5)) > ind(5))
# parallel trial design intervention indicator
df$int <- 0
df[df$cl > 3, 'int'] <- 1
# specify parameter values in the call for the data simulation below
des <- Model$new(
formula= ~ factor(t) + int - 1 +(1|gr(cl)*ar0(t)),
covariance = c(0.05,0.7),
mean = c(rep(0,5),0.2),
data = df,
family = gaussian()
)
ysim <- des$sim_data() # simulate some data from the model
fit1 <- des$MCML(y = ysim) # Default model fitting with SAEM
# use MCNR instead
fit2 <- des$MCML(y = ysim, method="mcnr")
# Non-linear model fitting example using the example provided by nlmer in lme4
data(Orange, package = "lme4")
# the lme4 example:
# startvec <- c(Asym = 200, xmid = 725, scal = 350)
# (nm1 <- lme4::nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
# Orange, start = startvec))
Orange <- as.data.frame(Orange)
Orange$Tree <- as.numeric(Orange$Tree)
# Here we can specify the model as a function.
model <- Model$new(
circumference ~ Asym/(1 + exp((xmid - (age))/scal)) - 1 + (Asym|gr(Tree)),
data = Orange,
family = gaussian(),
mean = c(200,725,350),
covariance = c(500),
var_par = 50
)
# for this example, we will use MCEM with adaptive MCMC sample sizes
nfit <- model$MCML(method = "mcem.adapt", iter.sampling = 1000)
summary(nfit)
}
Method fit()
MCML model fitting with the fastest options
Uses double Newton-Raphson method (with or without REML for Gaussian models). For details on the algorithm see MCML(). This function
uses the fastest set of options, including specialised model fitting for linear models. Note that no random effect
samples are drawn for Gaussian models, but can be subsequently drawn using mcmc_sample(). It is recommended to use the log
version of the covariance functions with this method as the Newton-Raphson steps can lead to negative values otherwise.
Usage
Model$fit(
niter = 100,
max_iter = 30,
tol = ifelse(self$family[[1]] == "gaussian" & self$family[[2]] == "identity", 1e-06,
10),
hist = 5,
k0 = 10,
reml = TRUE
)
Arguments
niterInteger. Number of samples for the random effects, ignored for Gaussian models, see examples.
max_iter
Integer. Maximum number of iterations.
tol
Scalar. The tolerance for the convergence criterion. For GLMMs this is the tolerance for
the Bayes Factor convergence criterion, for Gaussian linear models the tolerance is the difference
in the log-likelihood between successive iterations.
hist
Integer. The length of the running mean for the convergence criterion for non-Gaussian models.
k0
Integer. The expected number of iterations until convergence.
reml
Bool. For Gaussian models, whether to use REML or not.
Returns
A mcml model fit object
Examples
# Simulated trial data example using REML
# specify covariance starting values to potential prevent fit failure with random intialisation
data(SimTrial,package = "glmmrBase")
fit1 <- Model$new(
formula = y ~ int + factor(t) - 1 + (1|grlog(cl)*ar0log(t)),
data = SimTrial,
family = gaussian(),
covariance = log(c(0.05,0.8)),
)$fit(reml = TRUE)# Salamanders data example
data(Salamanders,package="glmmrBase")
# specify covariance starting values to prevent potential fit failure with random intialisation
model <- Model$new(
mating~fpop:mpop-1+(1|grlog(mnum))+(1|grlog(fnum)),
data = Salamanders,
family = binomial(),
covariance = c(-0.5, -0.5)
)
fit2 <- model$fit()
# Example using simulated data
#create example data with six clusters, five time periods, and five people per cluster-period
df <- nelder(~(cl(20)*t(10)) > ind(5))
# parallel trial design intervention indicator
df$int <- 0
df[df$cl > 10, 'int'] <- 1
# specify parameter values in the call for the data simulation below
des <- Model$new(
formula= ~ factor(t) + int - 1 +(1|grlog(cl)*ar0log(t)),
covariance = log(c(0.15,0.7)),
mean = c(rep(0,10),0.2),
data = df,
family = binomial()
)
ysim <- des$sim_data() # simulate some data from the model
des$update_y(ysim)
fit2 <- des$fit()
Method sparse()
Set whether to use sparse matrix methods for model calculations and fitting.
By default the model does not use sparse matrix methods.
Usage
Model$sparse(sparse = TRUE)
Arguments
sparseLogical indicating whether to use sparse matrix methods
Returns
None, called for effects
Method mcmc_sample()
Generate an MCMC sample of the random effects
Usage
Model$mcmc_sample(mcmc.pkg = "rstan", scaled = TRUE, constr.zero = 1)
Arguments
mcmc.pkgString. Either analytic for importance sampling from Gaussian posterior proposal, cmdstan for cmdstan
(requires the package cmdstanr), rstan to use rstan sampler (the default)
scaled
Logical. The random effects are sampled from an N(0,I) distribution. If TRUE this function returns the random effects rescaled to N(0,D), otherwise it returns the original samples.
constr.zero
Scalar. A Soft sum-to-zero constraint can be forced on the random effects so that their sum is N(0,constr.zero*Q). Small values, e.g. 0.001
may be useful if there is possible identifiability issues for intercept terms, such as in more complex, or higher dimensional, random effects structures like spatial models.
Returns
A matrix of samples of the random effects
Method importance_weights()
Returns the importance weights for the random effect samples. If using MCMC then weights are all equal.
Usage
Model$importance_weights()
Returns
A vector of the weights
Method gradient()
The gradient of the log-likelihood with respect to either the random effects or
the model parameters. The random effects are on the N(0,I) scale, i.e. scaled by the
Cholesky decomposition of the matrix D. To obtain the random effects from the last
model fit, see member function $u
Usage
Model$gradient(y, u, beta = FALSE)
Arguments
y(optional) Vector of outcome data, if not specified then data must have been set in another function.
u
(optional) Vector of random effects scaled by the Cholesky decomposition of D
beta
Logical. Whether the log gradient for the random effects (FALSE) or for the linear predictor parameters (TRUE)
Returns
A vector of the gradient
Method partial_sigma()
The partial derivatives of the covariance matrix Sigma with respect to the covariance
parameters. The function returns a list in order: Sigma, first order derivatives, second
order derivatives. The second order derivatives are ordered as the lower-triangular matrix
in column major order. Letting 'd(i)' mean the first-order partial derivative with respect
to parameter i, and d2(i,j) mean the second order derivative with respect to parameters i
and j, then if there were three covariance parameters the order of the output would be:
(sigma, d(1), d(2), d(3), d2(1,1), d2(1,2), d2(1,3), d2(2,2), d2(2,3), d2(3,3)).
Usage
Model$partial_sigma()
Returns
A list of matrices, see description for contents of the list.
Method u()
Returns the sample of random effects from the last model fit, or updates the samples for the model.
Usage
Model$u(scaled = TRUE, u)
Arguments
scaledLogical indicating whether the samples are on the N(0,I) scale (scaled=FALSE) or
N(0,D) scale (scaled=TRUE)
u
(optional) Matrix of random effect samples. If provided then the internal samples are replaced with these values. These samples should be N(0,I).
Returns
A matrix of random effect samples
Method log_likelihood()
The log likelihood for the GLMM. The random effects can be left
unspecified. If no random effects are provided, and there was a previous model fit with the same data y
then the random effects will be taken from that model. If there was no
previous model fit then the random effects are assumed to be all zero.
Usage
Model$log_likelihood(y, u)
Arguments
yA vector of outcome data
u
An optional matrix of random effect samples. This can be a single column.
Returns
The log-likelihood of the model parameters
Method calculator_instructions()
Prints the internal instructions and data used to calculate the linear predictor.
Internally the class uses a reverse polish notation to store and
calculate different functions, including user-specified non-linear mean functions. This
function will print all the steps. Mainly used for debugging and determining how the
class has interpreted non-linear model specifications.
Usage
Model$calculator_instructions()
Returns
None. Called for effects.
Method marginal()
Calculates the marginal effect of variable x. There are several options for
marginal effect and several types of conditioning or averaging. The type of marginal
effect can be the derivative of the mean with respect to x (dydx), the expected
difference E(y|x=a)-E(y|x=b) (diff), or the expected log ratio log(E(y|x=a)/E(y|x=b)) (ratio).
Other fixed effect variables can be set at specific values (at), set at their mean values
(atmeans), or averaged over (average). Averaging over a fixed effects variable here means
using all observed values of the variable in the relevant calculation.
The random effects can similarly be set at their
estimated value (re="estimated"), set to zero (re="zero"), set to a specific value
(re="at"), or averaged over (re="average"). Estimates of the expected values over the random
effects are generated using MCMC samples. MCMC samples are generated either through
MCML model fitting or using mcmc_sample. In the absence of samples average and estimated
will produce the same result. The standard errors are calculated using the delta method with one
of several options for the variance matrix of the fixed effect parameters.
Several of the arguments require the names
of the variables as given to the model object. Most variables are as specified in the formula,
factor variables are specified as the name of the variable_value, e.g. t_1. To see the names
of the stored parameters and data variables see the member function names().
Usage
Model$marginal(
x,
type,
re,
se,
at = c(),
atmeans = c(),
average = c(),
xvals = c(1, 0),
atvals = c(),
revals = c(),
oim = FALSE
)
Arguments
xString. Name of the variable to calculate the marginal effect for.
type
String. Either dydx for derivative, diff for difference, or ratio for log ratio. See description.
re
String. Either estimated to condition on estimated values, zero to set to zero, at to
provide specific values, or average to average over the random effects.
se
String. Type of standard error to use, either GLS for the GLS standard errors, KR for
Kenward-Roger estimated standard errors, KR2 for the improved Kenward-Roger correction (see small_sample_correction()),
or robust to use a robust sandwich estimator.
at
Optional. A vector of strings naming the fixed effects for which a specified value is given.
atmeans
Optional. A vector of strings naming the fixed effects that will be set at their mean value.
average
Optional. A vector of strings naming the fixed effects which will be averaged over.
xvals
Optional. A vector specifying the values of a and b for diff and ratio. The default is (1,0).
atvals
Optional. A vector specifying the values of fixed effects specified in at (in the same order).
revals
Optional. If re="at" then this argument provides a vector of values for the random effects.
oim
Logical. If TRUE use the observed information matrix, otherwise use the expected information matrix for standard error and related calculations.
Returns
A named vector with elements margin specifying the point estimate and se giving the standard error.
Method update_y()
Updates the outcome data y
Some functions require outcome data, which is by default set to all zero if no model fitting function
has been run. This function can update the interval y data.
Arguments
yVector of outcome data
Returns
None. Called for effects
Method set_trace()
Controls the information printed to the console for other functions.
Usage
Model$set_trace(trace)
Arguments
traceInteger, either 0 = no information, 1 = some information, 2 = all information
Returns
None. Called for effects.
Method clone()
The objects of this class are cloneable with this method.
Usage
Model$clone(deep = FALSE)
Arguments
deepWhether to make a deep clone.