geeDREstimation {CRTgeeDR} | R Documentation |

## Doubly Robust Inverse Probability Weighted Augmented GEE Estimator

### Description

This function implements a GEE estimator. It implements classical GEE, IPW-GEE, augmented GEE and IPW-Augmented GEE (Doubly robust).

### Usage

```
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)
```

### Arguments

`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 |

`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. |

### Details

The estimator is founds by solving:

` 0= \sum_{i=1}^M \Bigg[ \boldsymbol D_i^T \boldsymbol V_i^{-1} \boldsymbol W_i(\boldsymbol X_i, A_i, \boldsymbol \eta_W) \left( \boldsymbol Y_i - \boldsymbol B(\boldsymbol X_i, A_i, \boldsymbol \eta_B) \right) `

` \qquad + \sum_{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 \eta_B) -\boldsymbol \mu_i(\boldsymbol \beta,A_i=a)\Big) \Bigg]`

where `\boldsymbol D_i=\frac{\partial \boldsymbol \mu_i(\boldsymbol \beta,A_i)}{\partial \boldsymbol \beta^T}`

is the design matrix, `\boldsymbol V_i`

is the covariance matrix equal to `\boldsymbol U_i^{1/2} \boldsymbol C(\boldsymbol \alpha)\boldsymbol U_i^{1/2}`

with `\boldsymbol U_i`

a diagonal matrix with elements `{\rm var}(y_{ij})`

and `\boldsymbol C(\boldsymbol \alpha)`

is the working correlation structure with non-diagonal terms `\boldsymbol \alpha`

.
Parameters `\boldsymbol \alpha`

are estimated using simple moment estimators from the Pearson residuals.
The matrix of weights `\boldsymbol W_i(\boldsymbol X_i, A_i, \boldsymbol \eta_W)=diag\left[R_{ij}/\pi_{ij}(\boldsymbol X_i, A_i, \boldsymbol \eta_W)\right]_{j=1,\dots,n_{i}}`

, where `\pi_{ij}(\boldsymbol X_i, A_i, \boldsymbol \eta_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 \eta_B)`

, which is called the Outcome Model (OM), is a function linking `Y_{ij}`

with `\boldsymbol X_i`

and `A_i`

.
The `\boldsymbol \eta_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{\beta}_{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 \eta}_W`

in `\boldsymbol W_i(\boldsymbol X_i, A_i, \hat{ \boldsymbol \eta}_W)`

are obtained using a logistic regression and `\hat{ \boldsymbol \eta}_B`

in `\boldsymbol B(\boldsymbol X_i,A_i,\hat{ \boldsymbol \eta}_B)`

are obtained using a linear regression.

The variance of `\hat{\boldsymbol \beta}_{aug}`

is estimated by the sandwich variance estimator.
There are two external sources of variability that need to be accounted for: estimation of `\boldsymbol \eta_W`

for the PS and of `\boldsymbol \eta_B`

for the OM.
We denote `\boldsymbol \Omega=(\boldsymbol \beta, \boldsymbol \eta_W,\boldsymbol \eta_B)`

the estimated parameters of interest and nuisance parameters.
We can stack estimating functions and score functions for `\boldsymbol \Omega`

:

`\small \boldsymbol U_i(\boldsymbol \Omega)= \left( \begin{array}{c} \boldsymbol \Phi_i(\boldsymbol Y_i,\boldsymbol X_i,A_i,\boldsymbol \beta, \boldsymbol \eta_W, \boldsymbol \eta_B) \\ \boldsymbol S^W_i(\boldsymbol X_i, A_i, \boldsymbol \eta_W)\\ \boldsymbol S^B_i(\boldsymbol X_i, A_i, \boldsymbol \eta_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 \eta_W`

and `\boldsymbol \eta_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 \Omega)={{E\left[\frac{\partial \boldsymbol U_i(\boldsymbol \Omega)}{\partial \boldsymbol \Omega}\right]}^{-1}}^{T} \underbrace{{E\left[ \boldsymbol U_i(\boldsymbol \Omega)\boldsymbol U_i^T(\boldsymbol \Omega) \right]}}_{\boldsymbol \Delta_{adj}} \underbrace{E\left[\frac{\partial \boldsymbol U_i(\boldsymbol \Omega)}{\partial \boldsymbol \Omega}\right]^{-1} }_{\boldsymbol \Gamma^{-1}_{adj}}.`

### Value

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

### Author(s)

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.]

### References

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.

### Examples

```
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)
```

*CRTgeeDR*version 2.0.1 Index]