A model for case-control studies with optional prior distributions for the
coefficients, intercept, and auxiliary parameters.
stan_clogit( formula, data, subset, na.action = NULL, ..., strata, prior = normal(autoscale = TRUE), prior_covariance = decov(), prior_PD = FALSE, algorithm = c("sampling", "optimizing", "meanfield", "fullrank"), adapt_delta = NULL, QR = FALSE, sparse = FALSE )
Arguments
| formula, data, subset, na.action | Same as for |
||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| ... | Further arguments passed to the function in the rstan
package ( |
||||||||||
| strata | A factor indicating the groups in the data where the number of
successes (possibly one) is fixed by the research design. It may be useful
to use |
||||||||||
| prior | The prior distribution for the (non-hierarchical) regression coefficients. The default priors are described in the vignette
Prior
Distributions for rstanarm Models.
If not using the default,
See the priors help page for details on the families and
how to specify the arguments for all of the functions in the table above.
To omit a prior ---i.e., to use a flat (improper) uniform prior---
Note: Unless |
||||||||||
| prior_covariance | Cannot be |
||||||||||
| prior_PD | A logical scalar (defaulting to |
||||||||||
| algorithm | A string (possibly abbreviated) indicating the
estimation approach to use. Can be |
||||||||||
| adapt_delta | Only relevant if |
||||||||||
| QR | A logical scalar defaulting to |
||||||||||
| sparse | A logical scalar (defaulting to |
Value
A stanreg object is returned
for stan_clogit.
Details
The stan_clogit function is mostly similar in syntax to
clogit but rather than performing maximum
likelihood estimation of generalized linear models, full Bayesian
estimation is performed (if algorithm is "sampling") via
MCMC. The Bayesian model adds priors (independent by default) on the
coefficients of the GLM.
The data.frame passed to the data argument must be sorted by
the variable passed to the strata argument.
The formula may have group-specific terms like in
stan_glmer but should not allow the intercept to vary by the
stratifying variable, since there is no information in the data with which
to estimate such deviations in the intercept.
See also
stanreg-methods and
clogit.
The vignette for Bernoulli and binomial models, which has more
details on using stan_clogit.
http://mc-stan.org/rstanarm/articles/
Examples
dat <- infert[order(infert$stratum), ] # order by strata post <- stan_clogit(case ~ spontaneous + induced + (1 | education), strata = stratum, data = dat, subset = parity <= 2, QR = TRUE, chains = 2, iter = 500) # for speed only#> #> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 9.7e-05 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.97 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 500 [ 0%] (Warmup) #> Chain 1: Iteration: 50 / 500 [ 10%] (Warmup) #> Chain 1: Iteration: 100 / 500 [ 20%] (Warmup) #> Chain 1: Iteration: 150 / 500 [ 30%] (Warmup) #> Chain 1: Iteration: 200 / 500 [ 40%] (Warmup) #> Chain 1: Iteration: 250 / 500 [ 50%] (Warmup) #> Chain 1: Iteration: 251 / 500 [ 50%] (Sampling) #> Chain 1: Iteration: 300 / 500 [ 60%] (Sampling) #> Chain 1: Iteration: 350 / 500 [ 70%] (Sampling) #> Chain 1: Iteration: 400 / 500 [ 80%] (Sampling) #> Chain 1: Iteration: 450 / 500 [ 90%] (Sampling) #> Chain 1: Iteration: 500 / 500 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 0.282441 seconds (Warm-up) #> Chain 1: 0.088467 seconds (Sampling) #> Chain 1: 0.370908 seconds (Total) #> Chain 1: #> #> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2). #> Chain 2: #> Chain 2: Gradient evaluation took 4.4e-05 seconds #> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds. #> Chain 2: Adjust your expectations accordingly! #> Chain 2: #> Chain 2: #> Chain 2: Iteration: 1 / 500 [ 0%] (Warmup) #> Chain 2: Iteration: 50 / 500 [ 10%] (Warmup) #> Chain 2: Iteration: 100 / 500 [ 20%] (Warmup) #> Chain 2: Iteration: 150 / 500 [ 30%] (Warmup) #> Chain 2: Iteration: 200 / 500 [ 40%] (Warmup) #> Chain 2: Iteration: 250 / 500 [ 50%] (Warmup) #> Chain 2: Iteration: 251 / 500 [ 50%] (Sampling) #> Chain 2: Iteration: 300 / 500 [ 60%] (Sampling) #> Chain 2: Iteration: 350 / 500 [ 70%] (Sampling) #> Chain 2: Iteration: 400 / 500 [ 80%] (Sampling) #> Chain 2: Iteration: 450 / 500 [ 90%] (Sampling) #> Chain 2: Iteration: 500 / 500 [100%] (Sampling) #> Chain 2: #> Chain 2: Elapsed Time: 0.336151 seconds (Warm-up) #> Chain 2: 0.094005 seconds (Sampling) #> Chain 2: 0.430156 seconds (Total) #> Chain 2:#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#bulk-essnd <- dat[dat$parity > 2, c("case", "spontaneous", "induced", "education", "stratum")] # next line would fail without case and stratum variables pr <- posterior_epred(post, newdata = nd) # get predicted probabilities # not a random variable b/c probabilities add to 1 within strata all.equal(rep(sum(nd$case), nrow(pr)), rowSums(pr))#> [1] TRUE