geelm.RdEstimate mean structure parameters and their corresponding standard errors for generalized linear models with clustered or correlated observations by use of generalized estimating equations.
geelm.fit(x, y, id, offset, family, weights, control, corstr, start = NULL)
geelm(
formula,
id = NULL,
waves = NULL,
data = parent.frame(),
family = gaussian,
corstr = "independence",
Mv = 1,
weights = NULL,
corr.mat = NULL,
offset = NULL,
engine = "geeasy",
output = "geelm",
control = geelm.control()
)For glm:
logical values indicating whether the response vector and model
matrix used in the fitting process should be returned as components
of the returned value.
For glm.fit: x is a design matrix of dimension
n * p, and y is a vector of observations of length
n.
A vector identifying the clusters. If NULL, then each observation is assigned its own cluster.
this can be used to specify an a priori known
component to be included in the linear predictor during fitting.
This should be NULL or a numeric vector of length equal to
the number of cases. One or more offset terms can be
included in the formula instead or as well, and if more than one is
specified their sum is used. See model.offset.
A description of the error distribution and link function to be used
in the model. The argument can be one of three options: a family object,
a character string, or a list of functions. For more information on how to use family
objects, see Details below.
an optional vector of ‘prior weights’ to be used
in the fitting process. Should be NULL or a numeric vector.
A list of parameters for controlling the fitting process.
A character string specifying the correlation structure.
The default is "independence". Allowed structures are: "independence",
"exchangeable", "ar1", "m-dependent", "unstructured",
"fixed", and "userdefined". Any unique substring may be supplied.
If "fixed" or "userdefined", then corr.mat must be
specified. If "m-dependent", then Mv is relevant.
starting values for the parameters in the linear predictor.
A formula expression similar to that for glm,
An numeric vector identifying the time ordering within clusters
(i.e. levels of id). By default, data are assumed
to be sorted such that observations in a cluster are in consecutive rows
and higher numbered rows in a cluster are assumed to be later. Note that only the
ordering of the values in waves is used, NOT the numeric values themselves.
This means that e.g. having waves equal to c(1, 2, 3)
or c(1, 2, 7) within a cluster results in the same model.
An optional data frame containing the variables in the model.
For "m-dependent", the value for m.
The correlation matrix for "fixed". Matrix should
be symmetric with dimensions >= the maximum cluster size. If the correlation
structure is "userdefined", then this is a matrix describing which
correlations are the same.
Engine used to fit the model. The default, "geeasy" uses this
package (built on the geeM package), while "geepack" uses
the function geeglm from geepack to fit the model. Note that if
the geepack engine is used, the data are sorted according to id (and possibly
waves within id) and NAs are dropped before the data is used
(this differs from the standard in geepack).
Output object type. There are two options; 1) "geelm" (default), resulting in
an output that inherits the structure of geepacks geeglm object, or 2)
"geem" (or its alias "geeM") which results in an output that has the structure
of geeMs geem object.
An object of class geelm (inherits from geeglm) representing the fit.
It contains the following slots:
$coefficients: Coefficients from the mean structure model (betas) on their
original scales
$residuals: Pearson residuals, in the order of the inputted dataset (with NAs omitted).
$fitted.values: Fitted values (response scale), in the order of the inputted dataset
(with NAs omitted).
$rank: The rank of the model matrix, i.e. the number of estimated mean structure
coefficients.
$qr: QR decomposition of the model matrix (NA omitted).
$family: A family object specifying which exponential family was used for fitting
the mean structure model, see family for more information.
$linear.predictors: The linear predictor on the original scale.
$weights: Weights used for computations, in the order of the inputted dataset
(NAs omitted).
$prior.weights: The original weights used to produce this geeglm object (set
by user or defaulted to 1 for all observations).
$df.residuals: Residual degrees of freedom.
$y: Outcome variable, in the order of the inputted dataset (NAs omitted).
$model: The model.frame, ordered as the original inputted data with NAs omitted.
$call: The original function call that produced this geeglm object.
$formula: The formula used in the original call.
$terms: The terms of the formula used in the original call.
$data: The original dataset that was used for producing this geeglm object.
$offset: Offset used for fitting the model, ordered as the
original inputted data with NAs omitted.
$control: Value of control parameters used for fitting the
model.
$method: Internal function used for fitting the model.
$contrasts: Contrasts used in the model matrix.
$xlevels: Levels of factor variables used in the model
formula (if any).
$geese: An object containing further information about the
variance estimation, including a variance matrix for the
beta-coefficients ($vbeta), the estimated coefficients for
the working correlation matrix ($alpha), the estimated
dispersion parameter ($gamma), and the individual cluster
sizes ($clusz). See geese for more
information.
$modelInfo: Information about the link functions used for
fitting the mean, variance and scale structures of the model.
$id: IDs used for identifying the clusters, ordered as the
original inputted data with NAs omitted.
$corstr: Name of the correlation structured imposed on the
model. If the correlation structure requires further information,
it is stored in a suitably named attribute. For example, for
m-dependent correlation structures, the m scalar is available in an
attribute named Mv.
$cor.link: Link function used for the correlation structure.
$std.err: Method used to estimate the standard error of the mean structure
coefficients (betas).
Users may specify functions for link and variance functions, but the functions must be vectorized functions.
Offsets can be specified in the model formula, as in glm() or they may be
specified using the offset argument. If offsets are specified in both ways,
their sum is used as an offset.
For the "userdefined" correlation option, the function accepts a
matrix with consecutive integers. Each such integer represent a distinct
parameter that will be estimated. All entries given as 1 will be assumed
to be the same as each other and will be assumed to be possibly different
from entries with a 2, and so on.geelm only looks at the upper
triangle of the matrix. Any entry given as 0 will be fixed at 0.
If observations are dropped because they have a weight of 0, then the
denominator for the moment estimates of the correlation matrices are
calculated using the number of non-zero Pearson residuals for the
correlation structures unstructured, userdefined and
m-dependent with Mv>1. Therefore, residuals numerically
equal to 0 may cause problems in the calculation of correlation parameters.
Concerning the family argument: If the supplied argument is a character
string, then the string should correspond to one of the family objects.
In order to define a link function, a list must be created with the
components (LinkFun, VarFun, InvLink, InvLinkDeriv), all of which are
vectorized functions. If the components in the list are not named
as (LinkFun, VarFun, InvLink, InvLinkDeriv), then geelm
assumes that the functions are given in that order. LinkFun and VarFun
are the link and variance functions. InvLink and InvLinkDeriv are the inverse
of the link function and the derivative of the inverse of the link function
and so are decided by the choice of the link function.
geelm.fit():
# load data
data("respiratory")
respiratory$useid <- interaction(respiratory$center, respiratory$id)
# fit model
m <- geelm(outcome ~ treat + sex + age + baseline,
data = respiratory, id = useid,
family = "binomial", corstr = "exchangeable")
if (FALSE) { # \dontrun{
get_jack_se <- function(object, dat){
parm <- sapply(1:nrow(dat),
function(i){
dat.i <- dat[-i,]
coef(update(object, data=dat.i))
})
parm <- t(parm)
parm.mean <- apply(parm, 2, mean)
parm.cent <- sapply(1:nrow(parm),
function(i){
parm[i, ] - parm.mean
})
parm.cent <- t(parm.cent)
jack.var <- ((nrow(dat)-1) / nrow(dat)) * t(parm.cent) %*% parm.cent
jack.se <- sqrt(diag(jack.var))
jack.se
}
# load data
data("respiratory")
respiratory$useid <- interaction(respiratory$center, respiratory$id)
# fit model
obj <- geelm(outcome ~ treat + sex + age + baseline,
data = respiratory, id = useid,
family = "binomial", corstr = "exchangeable")
dat <- respiratory
get_jack_se(obj, dat)
summary(obj) |> coef()
} # }