View source: R/CRTgeeDR-main.R
geeDREstimation | R Documentation |
This function implements a GEE estimator. It implements classical GEE, IPW-GEE, augmented GEE and IPW-Augmented GEE (Doubly robust).
geeDREstimation(formula, id, data = parent.frame(), family = gaussian, corstr = "independence", Mv = 1, weights = NULL, aug = NULL, pi.a = 1/2, corr.mat = NULL, init.beta = NULL, init.alpha = NULL, init.phi = 1, scale.fix = FALSE, sandwich = TRUE, maxit = 20, tol = 1e-05, print.log = FALSE, typeweights = "VW", nameTRT = "TRT", model.weights = NULL, model.augmentation.trt = NULL, model.augmentation.ctrl = NULL, stepwise.augmentation = FALSE, stepwise.weights = FALSE, nameMISS = "MISSING", nameY = "OUTCOME", sandwich.nuisance = FALSE, fay.adjustment = FALSE, fay.bound = 0.75)
formula |
an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
id |
a vector which identifies the clusters. The length of "id" should be the same as the number of observations. Data are assumed to be sorted so that observations on a cluster are contiguous rows for all entities in the formula. |
data |
an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which CRTgeeDR is called. |
family |
a description of the error distribution and link function to be used in the model. This can be a character string naming a family function, a family function or the result of a call to a family function. (See family for details of family functions.) |
corstr |
a character string specifying the correlation structure. The following are permitted: '"independence"', '"exchangeable"', '"ar1"', '"unstructured"' and '"userdefined"' |
Mv |
for "m-dependent", the value for m |
weights |
A vector of weights for each observation. If an observation has weight 0, it is excluded from the calculations of any parameters. Observations with a NA anywhere (even in variables not included in the model) will be assigned a weight of 0. |
aug |
A list of vector (one for A=1 treated, one for A=0 control) for each observation representing E(Y|X,A=a). |
pi.a |
A number, Probability of treatment attribution P(A=1) |
corr.mat |
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. |
init.beta |
an optional vector with the initial values of beta. If not specified, then the intercept will be set to InvLink(mean(response)). init.beta must be specified if not using an intercept. |
init.alpha |
an optional scalar or vector giving the initial values for the correlation. If provided along with Mv>1 or unstructured correlation, then the user must ensure that the vector is of the appropriate length. |
init.phi |
an optional initial overdispersion parameter. If not supplied, initialized to 1. |
scale.fix |
if set to TRUE, then the scale parameter is fixed at the value of init.phi. |
sandwich |
if set to TRUE, the sandwich variance is provided together with the naive estimator of variance. |
maxit |
maximum number of iterations. |
tol |
tolerance in calculation of coefficients. |
print.log |
if set to TRUE, a report is printed. |
typeweights |
a character string specifying the weights implementation. The following are permitted: "GENMOD" for W^{1/2}V^{-1}W^{1/2}, "WV" for V^{-1}W |
nameTRT |
Name of the variable containing information for the treatment |
model.weights |
an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted for the propensity score. Must model the probability of being observed. |
model.augmentation.trt |
an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted for the ouctome model for the treated group (A=1). |
model.augmentation.ctrl |
an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted for the ouctome model for the control group (A=0). |
stepwise.augmentation |
if set to TRUE, a stepwise for the augmentation model is performed during the fit of the augmentation model for the OM |
stepwise.weights |
if set to TRUE, a stepwise for the propensity score is performed during the fit of the augmentation model for the OM |
nameMISS |
Name of the variable containing information for the Missing indicator |
nameY |
Name of the variable containing information for the outcome |
sandwich.nuisance |
if set to TRUE, the nuisance adjusted sandwich variance is provided. |
fay.adjustment |
if set to TRUE, the small-sample nuisance adjusted sandwich variance with Fay's adjustement is provided. |
fay.bound |
if set to 0.75 by default, bound value used for Fay's adjustement. |
The estimator is founds by solving:
0= ∑_{i=1}^M \Bigg[ \boldsymbol D_i^T \boldsymbol V_i^{-1} \boldsymbol W_i(\boldsymbol X_i, A_i, \boldsymbol η_W) ≤ft( \boldsymbol Y_i - \boldsymbol B(\boldsymbol X_i, A_i, \boldsymbol η_B) \right)
\qquad + ∑_{a=0,1} p^a(1-p)^{1-a} \boldsymbol D_i^T \boldsymbol V_i^{-1} \Big( \boldsymbol B(\boldsymbol X_i,A_i=a, \boldsymbol η_B) -\boldsymbol μ_i(\boldsymbol β,A_i=a)\Big) \Bigg]
where \boldsymbol D_i=\frac{\partial \boldsymbol μ_i(\boldsymbol β,A_i)}{\partial \boldsymbol β^T} is the design matrix, \boldsymbol V_i is the covariance matrix equal to \boldsymbol U_i^{1/2} \boldsymbol C(\boldsymbol α)\boldsymbol U_i^{1/2} with \boldsymbol U_i a diagonal matrix with elements {\rm var}(y_{ij}) and \boldsymbol C(\boldsymbol α) is the working correlation structure with non-diagonal terms \boldsymbol α. Parameters \boldsymbol α are estimated using simple moment estimators from the Pearson residuals. The matrix of weights \boldsymbol W_i(\boldsymbol X_i, A_i, \boldsymbol η_W)=diag≤ft[R_{ij}/π_{ij}(\boldsymbol X_i, A_i, \boldsymbol η_W)\right]_{j=1,…,n_{i}}, where π_{ij}(\boldsymbol X_i, A_i, \boldsymbol η_W)=P(R_{ij}|\boldsymbol X_i, A_i) is the Propensity score (PS). The function \boldsymbol B(\boldsymbol X_i,A_i=a,\boldsymbol η_B), which is called the Outcome Model (OM), is a function linking Y_{ij} with \boldsymbol X_i and A_i. The \boldsymbol η_B are nuisance parameters that are estimated. The estimator is most efficient if the OM is equal to E(\boldsymbol Y_i|\boldsymbol X_i,A_i=a) The estimator denoted \hat{β}_{aug} is found by solving the estimating equation. Although analytic solutions sometimes exist, coefficient estimates are generally obtained using an iterative procedure such as the Newton-Raphson method. Automatic implementation is such that, \hat{ \boldsymbol η}_W in \boldsymbol W_i(\boldsymbol X_i, A_i, \hat{ \boldsymbol η}_W) are obtained using a logistic regression and \hat{ \boldsymbol η}_B in \boldsymbol B(\boldsymbol X_i,A_i,\hat{ \boldsymbol η}_B) are obtained using a linear regression.
The variance of \hat{\boldsymbol β}_{aug} is estimated by the sandwich variance estimator. There are two external sources of variability that need to be accounted for: estimation of \boldsymbol η_W for the PS and of \boldsymbol η_B for the OM. We denote \boldsymbol Ω=(\boldsymbol β, \boldsymbol η_W,\boldsymbol η_B) the estimated parameters of interest and nuisance parameters. We can stack estimating functions and score functions for \boldsymbol Ω:
\small \boldsymbol U_i(\boldsymbol Ω)= ≤ft( \begin{array}{c} \boldsymbol Φ_i(\boldsymbol Y_i,\boldsymbol X_i,A_i,\boldsymbol β, \boldsymbol η_W, \boldsymbol η_B) \\ \boldsymbol S^W_i(\boldsymbol X_i, A_i, \boldsymbol η_W)\\ \boldsymbol S^B_i(\boldsymbol X_i, A_i, \boldsymbol η_B)\\ \end{array} \right)
where \boldsymbol S^W_i and \boldsymbol S^B_i represent the score equations for patients in cluster i for the estimation of \boldsymbol η_W and \boldsymbol η_B in the PS and the OM. A standard Taylor expansion paired with Slutzky's theorem and the central limit theorem give the sandwich estimator adjusted for nuisance parameters estimation in the OM and PS:
Var(\boldsymbol Ω)={{E≤ft[\frac{\partial \boldsymbol U_i(\boldsymbol Ω)}{\partial \boldsymbol Ω}\right]}^{-1}}^{T} \underbrace{{E≤ft[ \boldsymbol U_i(\boldsymbol Ω)\boldsymbol U_i^T(\boldsymbol Ω) \right]}}_{\boldsymbol Δ_{adj}} \underbrace{E≤ft[\frac{\partial \boldsymbol U_i(\boldsymbol Ω)}{\partial \boldsymbol Ω}\right]^{-1} }_{\boldsymbol Γ^{-1}_{adj}}.
An object of type 'CRTgeeDR'
$beta Final values for regressors estimates
$phi scale parameter estimate
$alpha Final values for association parameters in the working correlation structure when exchangeable
$coefnames Name of the regressors in the main regression
$niter Number of iteration done by the algorithm before convergence
$converged convergence status
$var.naiv Variance of the estimates model based (naive)
$var Variance of the estimates sandwich
$var.nuisance Variance of the estimates nuisance adjusted sandwich
$var.fay Variance of the estimates nuisance adjusted sandwich with Fay correction for small samples
$call Call function
$corr Correlation structure used
$clusz Number of unit in each cluster
$FunList List of function associated with the family
$X design matrix for the main regression
$offset Offset specified in the regression
$eta predicted values
$weights Weights vector used in the diagonal term for the IPW
$ps.model Summary of the regression fitted for the PS if computed internally
$om.model.trt Summary of the regression fitted for the OM for treated if computed internally
$om.model.ctrl Summary of the regression fitted for the OM for control if computed internally
Melanie Prague [based on R packages 'geeM' L. S. McDaniel, N. C. Henderson, and P. J. Rathouz. Fast Pure R Implementation of GEE: Application of the Matrix Package. The R Journal, 5(1):181-188, June 2013.]
Details regarding implementation can be found in
'Augmented GEE for improving efficiency and validity of estimation in cluster randomized trials by leveraging cluster-and individual-level covariates' - 2012 - Stephens A., Tchetgen Tchetgen E. and De Gruttola V. : Stat Med 31(10) - 915-930.
'Accounting for interactions and complex inter-subject dependency for estimating treatment effect in cluster randomized trials with missing at random outcomes' - 2015 - Prague M., Wang R., Stephens A., Tchetgen Tchetgen E. and De Gruttola V. : in revision.
'Fast Pure R Implementation of GEE: Application of the Matrix Package' - 2013 - McDaniel, Lee S and Henderson, Nicholas C and Rathouz, Paul J : The R Journal 5(1) - 181-197.
'Small-Sample Adjustments for Wald-Type Tests Using Sandwich Estimators' - 2001 - Fay, Michael P and Graubard, Barry I : Biometrics 57(4) - 1198-1206.
data(data.sim) ## Not run: #### STANDARD GEE geeresults<-geeDREstimation(formula=OUTCOME~TRT, id="CLUSTER" , data = data.sim, family = "binomial", corstr = "independence") summary(geeresults) #### IPW GEE ipwresults<-geeDREstimation(formula=OUTCOME~TRT, id="CLUSTER" , data = data.sim, family = "binomial", corstr = "independence", model.weights=I(MISSING==0)~TRT*AGE) summary(ipwresults) #### AUGMENTED GEE augresults<-geeDREstimation(formula=OUTCOME~TRT, id="CLUSTER" , data = data.sim, family = "binomial", corstr = "independence", model.augmentation.trt=OUTCOME~AGE, model.augmentation.ctrl=OUTCOME~AGE, stepwise.augmentation=FALSE) summary(augresults) ## End(Not run) #### DOUBLY ROBUST drresults<-geeDREstimation(formula=OUTCOME~TRT, id="CLUSTER" , data = data.sim, family = "binomial", corstr = "independence", model.weights=I(MISSING==0)~TRT*AGE, model.augmentation.trt=OUTCOME~AGE, model.augmentation.ctrl=OUTCOME~AGE, stepwise.augmentation=FALSE) summary(drresults)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.