*! version 1.0.0 05Jul2021

Title

stipw -- Inverse probability weighted parametric survival models with variance obtained via M-estimation

Syntax

stipw
(logit tvar tmvarlist
[, tmoptions])
[if] [in]
,
distribution(distname)
[options]

tvar must be a binary variable with 1 = treatment/exposure and 0 = control.

tmvarlist specifies the variables that predict treatment assignment/exposure in the treatment/exposure model.

tmoptions

Treatment/Exposure Model
noconstant suppress constant from treatment/exposure model

offset(varname) include varname in model with coefficient constrained to 1

tcoef displays coefficient table from treatment/exposure model

Maximization
maximize_options control the maximization process; seldom used

Outcome Model: streg
distribution(exponential) (weighted) exponential survival outcome model

distribution(weibull) (weighted) Weibull survival outcome model

distribution(gompertz) (weighted) Gompertz survival outcome model

distribution(loglogistic) (weighted) loglogistic survival outcome model

distribution(lognormal) (weighted) lognormal survival outcome model

ancillary specifies that tvar should be used to model ancillary parameter

ocoef displays coefficient table from outcome model, before variance is updated with M-estimation

oheader displays header from outcome model, before variance is updated with M-estimation

Outcome Model: stpm2
distribution(rp) (weighted) flexible parametric (Royston-Parmar) survival outcome model

bknots(knotslist) boundary knots for baseline

bknotstvc(knotslist) boundary knots for time-dependent treatment/exposure

df(#) degrees of freedom for baseline hazard function

dftvc(#) degrees of freedom for time-dependent treatment/exposure

failconvlininit automatically try lininit option if convergence fails

knots(numlist) knot locations for baseline hazard

knotstvc(numlist) knot locations for time-dependent treatment/exposure

knscale(scale) scale for user-defined knots (default scale is time)

noorthog do not use orthogonal transformation of splines variables

ocoef displays coefficient table from outcome model, before variance is updated with M-estimation

oheader displays header from outcome model, before variance is updated with M-estimation

SE/Robust
vce(vcetype) vcetype may be mestimation, the default, or robust

Advanced
ipwtype(string) type of IPW weight: stabilised (default) or unstabilsied

genweight(newvar) name of generated weight variable

genflag(newvar) name of generated flag variable that indicates which observations have been used in the analysis

stsetupdate specifies that the data in memory should be updated with stset that specifies the weights

Reporting
level(#) set confidence level; default is level(95)

noheader suppress header from coefficient table

nohr * do not report hazard ratios

tratio * report time ratios

noshow * do not show st setting information

eform ** exponentiate coefficients

alleq ** report all equations

keepcons ** do not drop constraints in ml routine (these can be seen using constraint dir)

showcons ** list constraints in output

display_options control columns and column formats and line width

Maximisation
lininit ** obtain initial values by first fitting a linear function of ln(time)

maximize_options control the maximization process; seldom used

coeflegend display legend instead of statistics

selegend display legend of standard errors instead of statistics

You must stset your data without weights and without the id option before using stipw; see
[ST] stset.

tmvarlist may not contain factor variables. You must create dummy variables when fitting the models.

* This option is only available for streg outcome models.

** This option is only available for stpm2 outcome models.

Description

stipw performs an inverse probability weighted (IPW) analysis on survival
data. The command fits a logistic regression model to model treatment/exposure (tvar against
confounders tmvarlist). By default, stabilised weights are calculated and therefore a second logistic
regression model is fit with no confounders. This is used to obtain the numerator for the weights (the treatment/exposure prevalence).
The second logistic regression model is not required for unstabilised weights.
The weights are named _stipw_weight, unless genweight(newvar) is specified.
The data is then stset with the weights (the original stset is preserved, unless stsetupdate is specified).
The specified distname survival model is fit to the weighted data. The program then estimates the variance
using a closed-form variance estimator obtained via M-estimation. This is the default variance estimator and the stored variance matrix is updated.
A robust variance estimator can be used instead, however, robust standard errors have been shown to be conservative in some scenarios (Austin 2016) and do
not appropriately take into account that the weights themselves are estimates and not fixed.

The variable _stipw_weight is created to store the weights. This is overwritten in
subsequent runs of stipw.

Survival model fitting commands supported include [ST] streg
and stpm2 (with scale(hazard) only).
Postestimation commands follow stipw
in the same way as they do for streg and stpm2, see stipw_postestimation.
The following are not supported with stipw:
generalized gamma models,
frailty and shared-frailty models,
relative survival models,
cure models and
multiple-record-per-subject survival data.

stipw with distribution(rp) creates new variables _rcs* and _d_rcs*.
If dftvc or knotstvc is specified, stipw with distribution(rp)
also creates new variables _rcs_tvar* and _d_rcs_tvar* where tvar is the
treatment/exposure variable. If there is delayed entry, new variables _s0_rcs* (and if applicable
_s0_rcs_tvar*) are created. These variables are all dropped in subsequent runs of stipw with
distribution(rp).

The option stsetupdate changes the data in memory in that it respecifies the stset command
with the weights, conditioning on the estimation sample (_stipw_flag).

Any observations with a missing tvar value or any missing values in tmvarlist
are excluded from the analysis. Any observations with _st == 0 are also excluded
from the analysis. A note is displayed if observations are excluded with the number excluded.
_stipw_flag is created to indicate which observations were used in the analysis.
_stipw_flag is overwritten in subsequent runs of stipw.

stipw requires package dm79, which can be installed using:

. net install dm79, from(http://www.stata.com/stb/stb56)

stipw with distribution(rp) also requires the stpm2 and
rcsgen packages, installed in the standard manner. Please ensure that you have installed
stpm2, version 1.7.5 May2021 or later.

Options

Treatment/Exposure Model

noconstant, offset(varname); see [R] estimation options.

tcoef displays coeffiecent table from logit treatment/exposure model(s). Default is to omit table.

Maximisation

lininit (only applicable to stpm2 models). This obtains initial values
for stpm2 by fitting only the first spline basis function (i.e. a linear function of log survival time).
This option is seldom needed.

maximize_options:
difficult,
technique(maximize##algorithm_spec:algorithm_spec),
iterate(#),
nolog,
trace,
gradient,
showstep,
hessian,
showtolerance,
tolerance(#),
ltolerance(#),
nrtolerance(#), and
nonrtolerance; see [R] maximize. These options are seldom used.

Outcome Model: streg

distribution(distname) option must be specified. This specifies the distribution of
the weighted survival model. The options in this sub-section follow for: distribution(exponential),
distribution(weibull), distribution(gompertz), distribution(lognormal)
and distribution(loglogistic).

ancillary specifies that tvar is used to model the ancillary parameter. By default, the ancillary
parameter does not depend on tvar.

ocoef displays coeffiecent table from streg outcome model, before the variance
is updated with M-estimation. Default is to omit table. This is an intermediate step and does not
correspond to the display options for the final coefficient table.

oheader displays header from streg outcome model, before the variance
is updated with M-estimation. Default is to omit header. This is an intermediate step and does not
correspond to the display options for the final coefficient table.

Outcome Model: stpm2

distribution(distname) option must be specified. This specifies the distribution of
the weighted survival model. The options in this sub-section follow for distribution(rp).

bknots(knotslist) knotslist is a two-element numlist giving
the boundary knots. By default these are located at the minimum and maximum
of the uncensored survival times. They are specified on the scale defined
by knscale().

bknotstvc(knotslist) knotslist is a two-element numlist giving
the boundary knots for the time-dependent treatment/exposure tvar.
By default these are the same as for the bknots option.
They are specified on the scale defined by knscale().
dftvc(#) or knotstvc(# [# ...]) need to be specified in order to
specify bknotstvc(knotslist).

df(#) specifies the degrees of freedom for the restricted
cubic spline function used for the baseline function. # must be between
1 and 10, but usually a value between 1 and 4 is sufficient, with 3 being a common choice.
Exactly one of knots() and df() must be specified.
The knots are placed at the following centiles of the
distribution of the uncensored log survival times:

df knots Centile positions

1 0 (no knots)
2 1 50
3 2 33 67
4 3 25 50 75
5 4 20 40 60 80
6 5 17 33 50 67 83
7 6 14 29 43 57 71 86
8 7 12.5 25 37.5 50 62.5 75 87.5
9 8 11.1 22.2 33.3 44.4 55.6 66.7 77.8 88.9
10 9 10 20 30 40 50 60 70 80 90

Note that these are interior knots and there are also boundary knots
placed at the minimum and maximum of the distribution of uncensored survival
times.

dftvc(#) gives the degrees of freedom for the time-dependent treatment/exposure tvar.
The potential degrees of freedom are listed under the df() option.
With 1 degree of freedom a linear effect of log time is fitted.
The knotstvc() option is not applicable if the dftvc() option
is specified.

failconvlininit automatically tries the lininit option of the
model fails to converge.

knots(# [# ...]) specifies knot locations for the baseline distribution
function, as opposed to the default locations set by df(). Note that
the locations of the knots are placed on the scale defined by knscale().
However, the scale used by the restricted cubic spline function is always
log time.

knotstvc(# [# ...]) specifies the location of the interior knots for
the time-dependent treatment/exposure tvar, as opposed to the default locations set by
dftvc().

knscale(scale) sets the scale on which user-defined knots are specified.
knscale(time) denotes the original time scale, knscale(log) the
log time scale and knscale(centile) specifies that the knots
are taken to be centile positions in the distribution of the uncensored log
survival times. The default is knscale(time).

noorthog suppresses orthogonal transformation of spline variables.

ocoef displays coeffiecent table from stpm2 outcome model, before the variance
is updated with M-estimation. Default is to omit table. This is an intermediate step and does not
correspond to the display options for the final coefficient table.

oheader displays header from stpm2 outcome model, before the variance
is updated with M-estimation. Default is to omit header. This is an intermediate step and does not
correspond to the display options for the final coefficient table. This option needs to be
specified with ocoef, as the header cannot be displayed without the coefficient table.

SE/Robust

vce(vcetype) specifies the type of standard error reported. The default is mestimation.
This provides a closed-form variance estimator, which appropriately takes into account that the
weights are estimated and not fixed. robust can be specified, but this will not take
into account that the weights are estimated.

Advanced

ipwtype(string) specifies the type of IPW weight to be used: stabilised (default) or unstabilised.
Stabilised weights require the estimation of the treatment/exposure prevalence and therefore there is an additional row/column in the full variance matrix, e(V_full).

genweight(newvar) specifies the name of the generated weight variable. By default, weights are stored in _stipw_weight.

genflag(newvar) specifies the name of the generated flag variable to show which observations have been used in the analysis.
By default, the indicator is stored in _stipw_flag.

stsetupdate changes the data in memory in that it respecifies the stset command
with the weights, conditioning on the estimation sample (_stipw_flag).

Reporting

level(#) specifies the confidence level, as a percentage, for confidence
intervals. The default is level(95) or as set by set level.

noheader suppresses the output header.

nohr (only applicable to streg models) specifies that coefficients
rather than exponentiated coefficients be displayed, that is, that coefficients
rather than hazard ratios be displayed. This option affects only how coefficients
are displayed, not how they are estimated. This option is valid only for models
with a natural proportional-hazards parameterization: exponential, Weibull and Gompertz.
These models, by default, report hazards ratios (exponentiated coefficients).

tratio (only applicable to streg models) specifies that
exponentiated coefficients, which are interpreted as time ratios, be displayed.
tratio is appropriate only for the loglogistic and lognormal or for the exponential
and Weibull models when fit in the accelerated failure-time metric.

noshow (only applicable to streg models) prevents stipw
from showing the key st variables. This option is rarely used because most people
type stset, show or stset, noshow to set once and for all whether
they want to see these variables mentioned at the top of the output of every st
command; see [ST] stset. stpm2 models must use noshow.

eform (only applicable to stpm2 models) reports the exponentiated coefficents.
This gives the hazard ratio of the treatment/exposure tvar if it is not time-dependent, i.e.,
if dftvc(#) and knotstvc(# [# ...]) are not specified.

keepcons (only applicable to stpm2 models) prevents the constraints
imposed by stpm2 on the derivatives of the spline function when fitting
delayed entry models being dropped. This can be viewed using constraint dir.
By default, the constraints are dropped.

alleq (only applicable to stpm2 models) reports all equations used by ml.
The models are fitted by using various constraints for parameters associated with
the derivatives of the spline functions. These parameters are generally not of
interest and thus are not shown by default. In addition, an extra equation is used
when fitting delayed entry models, and again this is not shown by default.

showcons (only applicable to stpm2 models). The constraints used by
stpm2 for the derivatives of the spline function and when fitting delayed
entry models are not listed by default. Use of this option lists them in the output.

display_options:
noci,
nopvalues,
cformat(%fmt),
pformat(%fmt),
sformat(%fmt),
nolstretch; see [R] estimation options.

The following options are available with stipw:

coeflegend; see [R] estimation options.

selegend specifies that the legend of the standard errors of the coefficients and how to specify them in an expression be displayed
rather than displaying the statistics for the standard errors of the coefficients.

Methods

The formulas below define unstabilised (uw) and stabilised (sw) weights, where A denotes
the binary treatment/exposure variable tvar, X denotes the confounders tmvarlist and I represents
the indicator function:

uw = I(A=1)/P(A=1|X) + I(A=0)/P(A=0|X)

sw = I(A=1)*P(A=1)/P(A=1|X) + I(A=0)*P(A=0)/P(A=0|X)

Point estimates are calculated from the weighted data in the usual way by maximum likelihood estimation with streg and stpm2.
The variance (of all model parameteres - treatment/exposure model(s) and outcome model)
is estimated using a closed-form variance estimator obtained via M-estimation.
M-estimation is chosen as the default variance estimator, as robust standard errors have been found to be conservative in some cases by Austin 2016
and do not appropriately take into account that the weights are estimated and not fixed.
M-estimation first involves defining a set of estimating equations, similar to those in Williamson et al 2014, but applied
to survival outcome models. The variance is then estimated as described by Stefanski and Boos 2002.
See Hill et al (in draft) for full details and formulas on how the variance is estimated.

For further reading, see Hajage et al 2018 and Shu et al 2020 for closed form variance estimators for IPW Cox models.

Example 1: streg model

Load example dataset:

. webuse brcancer, clear

stset the data:

. stset rectime, f(censrec==1) scale(365.24)

Exposure of interest is hormon.

Fit an unweighted Weibull model to the data (ignoring confounders) using streg:

. streg hormon, distribution(weibull)

We wish to adjust for confounders x1, x2, x3, x5, x6 and x7.

Perform an IPW analysis using stipw. Fit a weighted Weibull model (weights from a logit model) and obtain standard errors from M-estimation (this is the default):

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(weibull)

Repeat, but obtain robust standard errors for comparison:

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(weibull) vce(robust)

Repeat, but show the output for each stage:

. stipw (logit hormon x1 x2 x3 x5 x6 x7, tcoef) , distribution(weibull) ocoef oheader

Repeat and store the stablised weights (the default) in sweight and indicator in flag:

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(weibull) genweight(sweight) genflag(flag)

Now use unstabilised weights. Store the unstabilised weights in uweight:

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(weibull) ipwtype(unstabilised) genweight(uweight)

Fit the model using hormon to model the ancillary parameter:

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(weibull) ancillary

Repeat, but now update the stset command:

. char list

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(weibull) ancillary stsetupdate

. char list

You will now need to stset the data without weights if you wish to run another stipw command.

Repeat, reporting time ratios and set CI level as 90%:

. stset rectime, f(censrec==1) scale(365.24)

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(weibull) ancillary tratio level(90)

Example 2: stpm2 model

Load example dataset:

. webuse brcancer, clear

stset the data:

. stset rectime, f(censrec==1) scale(365.24)

Exposure of interest is hormon.

Fit an unweighted Royston-Parmar (RP) model to the data (ignoring confounders) with 3 degrees of freedom using stpm2:

. stpm2 hormon, scale(hazard) df(3)

We wish to adjust for confounders x1, x2, x3, x5, x6 and x7.

Perform an IPW analysis using stipw. Fit a weighted RP model (weights from a logit model) and obtain standard errors from M-estimation (this is the default):

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(rp) df(3)

Repeat, but obtain robust standard errors for comparison:

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(rp) df(3) vce(robust)

Repeat, but show the output for each stage:

. stipw (logit hormon x1 x2 x3 x5 x6 x7, tcoef) , distribution(rp) df(3) ocoef oheader

By default, stabilised weights are used. Now use unstabilised weights and store in uweight:

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(rp) df(3) ipwtype(unstabilised) genweight(uweight)

Repeat, but as hormon is not time-dependent, when we use the eform option this will be the hazard ratio:

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(rp) df(3) eform

Fit an unweighted RP model to the data (ignoring confounders) with time-dependent hormon:

. stpm2 hormon, scale(hazard) df(3) tvc(hormon) dftvc(2)

Fit the corresponding weighted model with stipw and obtain M-estimation standard errors

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(rp) df(3) dftvc(2)

Repeat, displaying all the equations and all the constraints

. stipw (logit hormon x1 x2 x3 x5 x6 x7) , distribution(rp) df(3) dftvc(2) alleq showcons

Stored results

stipw stores the following in e(), most of which come from streg/stpm2.

* This object is only available for streg outcome models.

** This object is only available for stpm2 outcome models.

Scalars

e(N) number of observations

e(n0) number of observations in the control group

e(n1) number of observations in the treatment/exposure group

e(N_sub) * number of subjects

e(N_fail) * number of failures

e(del_entry) ** indicator for delayed entry

e(k) number of parameters

e(k_eq) number of equations in e(b)

e(k_eq_model) number of equations in overall model test

e(k_aux) * number of auxiliary parameters

e(k_dv) number of dependent variables

e(df_m) * model degrees of freedom (includes ancillary parameter if specified)

e(dfbase) ** degrees of freedom in baseline hazard

e(df_tvar) ** degrees of freedom for time-dependent treatment/exposure variable tvar

e(nxbterms) ** number of xb terms (_rcs*) including time-dependent tvar terms

e(ndxbterms) ** number of dxb terms (_d_rcs*) including time-dependent tvar terms

e(ll) log likelihood

e(ll_0) * log likelihood, constant-only model

e(dev) ** deviance = -2 * log likelihood

e(AIC) ** AIC

e(AIC) ** BIC

e(chi2) * chi-squared test statistic simultaneously testing all occurences of tvar (including ancillary/splines)

e(risk) * total time at risk

e(aux_p) * ancillary parameter (weibull)

e(gamma) * ancillary parameter (gompertz, loglogistic)

e(sigma) * ancillary parameter (lnormal)

e(p) * p-value for model test

e(rank) rank of e(V)

e(ic) number of iterations

e(rc) return code for streg/stpm2 outcome model

e(rc_logit) return code for logit treatment/exposure model

e(rc_logit2) return code for the second logit model (no confounders) if stabilised weights are used

e(converged) 1 if converged, 0 otherwise (streg/stpm2 outcome model)

e(converged_logit) 1 if converged, 0 otherwise (logit treatment/exposure model)

e(converged_logit2) 1 if converged, 0 otherwise (second logit model (no confounders) if stabilised weights are used)

Macros

e(cmd) model (streg) or stpm2

e(cmd2) * streg

e(cmd3) stipw

e(cmdline) command as typed

e(cmdline_logit) command as used for logit treatment/exposure model

e(cmdline_logit2) command as used for the second logit model (no confounders) if stabilised weights are used

e(cmdline_streg) * command as used for streg outcome model.

e(cmdline_stpm2) ** command as used for stpm2 outcome model.

e(dead) * _d

e(depvar) _t (streg) or _d _t (stpm2)

e(tvar) name of treatment/exposure variable tvar

e(offset_logit) name of offset variable used in the treatment/exposure logit model

e(varnames) ** treatment/exposure variable tvar

e(varlist) ** treatment/exposure variable tvar

e(tvc) ** time-dependent variable: treatment/exposure variable tvar or missing

e(drcsterms_base) ** name of the _d_rcs* baseline hazard splines

e(rcsterms_base) ** name of the _rcs* baseline hazard splines

e(drcsterms_tvar) ** name of the _d_rcs_tvar* time-dependent treatment/exposure tvar splines

e(rcsterms_tvar) ** name of the _rcs_tvar* time-dependent treatment/exposure tvar splines

e(ln_bhknots) ** log of all the knot locations for the baseline hazard

e(bhknots) ** interior knot locations for the baseline hazard on the time scale

e(ln_tvcknots_tvar) ** log of all the knot locations for time-dependent treatment/exposure tvar variable

e(tvcknots_tvar) ** interior knot locations for time-dependent treatment/exposure tvar variable on the time scale

e(boundary_knots) ** boundary knot locations for the baseline hazard

e(boundary_knots_tvar) ** boundary knot locations for time-dependent treatment/exposure tvar variable

e(scale) ** the scale on which the survival model is fitted: hazard

e(orthog) ** indicator of orthogonal transformation of spline variables

e(title) * title in estimation output

e(wtype) weight type: pweight only

e(ipwtype) IPW weight type: stabilised or unstabilised

e(wexp) weight expression

e(t0) * _t0

e(vce) vcetype specified in vce()

e(vcetype) title used to label Std. Err.

e(frm2) * hazard or time

e(chi2type) * Wald; type of model chi-squared test

e(stcurve) * stcurve

e(opt) type of optimization

e(which) max or min; whether optimizer is to perform
maximization or minimization

e(ml_method) type of ml method

e(user) name of likelihood-evaluator program

e(technique) maximization technique

e(properties) b V

e(predict) program used to implement predict

e(predict_sub) * predict subprogram

Matrices

e(b) coefficient vector

e(b_full) full coefficient vector including estimates from treatment/exposure model (produced with vce(mestimation))

e(ilog) iteration log (up to 20 iterations)

e(gradient) gradient vector

e(Cns) ** constraints matrix

e(R_bh) ** matrix used for orthogonal transformation of baseline hazard spline variables

e(R_tvar) ** matrix used for orthogonal transformation of time-dependent treatment/exposure tvar spline variables

e(V) variance-covariance matrix of the estimators

e(V_full) full variance-covariance matrix including estimates from treatment/exposure model (produced with vce(mestimation))

e(V_A) A matrix used to calculate vce(mestimation) standard errors

e(V_B) B matrix used to calculate vce(mestimation) standard errors

e(V_modelbased) model-based variance

e(V_robust) standard robust variance

Functions

e(sample) marks estimation sample

Authors

Micki Hill

Biostatistics Research Group

Department of Health Sciences

University of Leicester

E-mail: mh594@le.ac.uk

Paul C. Lambert

Biostatistics Research Group

Department of Health Sciences

University of Leicester

and

Department of Medical Epidemiology and Biostatistics

Karolinska Institutet

Michael J. Crowther

Department of Medical Epidemiology and Biostatistics

Karolinska Institutet

Please report any errors you may find.

References

Austin PC. Variance estimation when using inverse probability of treatment weighting (IPTW) with survival analysis. Statistics in Medicine. 2016;35(30):5642-55.

Hajage D, Chauvet G, Belin L, Lafourcade A, Tubach F, De Rycke Y. Closed‐form variance estimator for weighted propensity score estimators with survival outcome. Biometrical Journal. 2018;60(6):1151-63.

Shu D, Young JG, Toh S, Wang R. Variance estimation in inverse probability weighted Cox models. Biometrics. 2020;1-17.

Stefanski LA and Boos DD. The calculus of M-estimation. The American Statistician. 2002;56;29-38.

Williamson EJ, Forbes A, White IR. Variance reduction in randomised trials by inverse probability weighting using the propensity score. Statistics in Medicine. 2014;33(5):721-37.