elrm (version 1.2.2)

elrm: elrm: exact-like inference in logistic regression models

Description

elrm implements a modification of the Markov Chain Monte Carlo algorithm proposed by Forster et al. (2003) to approximate exact conditional inference for logistic regression models. The modifications can handle larger datasets than the original algorithm (Zamar 2006). Exact conditional inference is based on the distribution of the sufficient statistics for the parameters of interest given the sufficient statistics for the remaining nuisance parameters. Using model formula notation, users specify a logistic model and model terms of interest for exact inference.

Usage

elrm(formula, interest, r = 4, iter = 1000, dataset, burnIn = 0, alpha = 0.05)

Arguments

formula

a formula object that contains a symbolic description of the logistic regression model of interest in the usual R formula format. One exception is that the binomial response should be specified as success/trials, where success gives the number of successes and trials gives the number of binomial trials for each row of dataset.

interest

a formula object that contains a symbolic description of the model terms for which exact conditional inference is of interest.

r

a parameter of the MCMC algorithm that influences how the Markov chain moves around the state space. Small values of r cause the chain to take small, relatively frequent steps through the state space; larger values cause larger, less frequent steps. The value of r must be an even integer less than or equal to the length of the response vector. Typical values are 4, 6 or 8; default=4.

iter

an integer representing the number of Markov chain iterations to make (must be larger than or equal to 1000); default=1000.

dataset

a data.frame object where the data are stored.

burnIn

the burn-in period to use when conducting inference. Values of the Markov chain in the burn-in period are discarded; default=0.

alpha

determines the level used for confidence intervals; default=0.05.

Value

coeffs

a vector containing the parameter estimates.

coeffs.ci

a list containing (1-alpha)*100% confidence intervals for each parameter of interest.

p.values

a vector containing the estimated p-value for jointly testing that the parameters of interest are simultaneously equal to zero, and the full conditional p-values from separately testing each parameter equal to zero.

p.values.se

a vector containing the Monte Carlo standard errors of the estimated p-values of each term of interest.

mc

an mcmc object containing the Markov chain of sampled values of the sufficient statistics for the parameters of interest. Columns correspond to parameters; rows to Monte Carlo iterations.

mc.size

a vector containing the lengths of the extracted Markov chains used in testing each parameter. The length of the Markov chain used for the joint test (i.e., iter) is also included as the first element.

obs.suff.stat

a vector containing the observed value of the sufficient statistic for each parameter of interest.

distribution

a list containing distribution tables for the sampled values of the sufficient statistic of the parameters of interest conditional on all the rest.

call.history

a list composed of the matched call and the history of calls to update().

dataset

the data.frame object that was passed to elrm() as an argument.

last

the last response vector sampled by the Markov chain.

r

the value of r passed to elrm() as an argument.

ci.level

the level used when constructing the confidence intervals for the parameters of interest. The level is calculated as (1-alpha)*100%.

Warning

The labels of the terms in the in the interest model should match those found in the formula model. Thus, the term.labels attribute of terms.formula(interest) should match those found in terms.formula(formula). Please see the Examples section for more details.

Details

The function summary() (i.e., summary.elrm) can be used to obtain or print a summary of the results. Each estimated exact p-value is based on the conditional probabilities test. The Monte Carlo standard error of each p-value is computed by the batch-means method (Geyer C.J. 1992). Inference on each parameter must be based on a Markov chain of at least 1000 iterations, otherwise NA is returned. If the observed value of the sufficient statistic for a parameter is either the maximum or the minimum value sampled, the MUE of the parameter is given instead of the CMLE. In such cases, the resulting confidence interval is open-ended on one side.

Apart from the documentation files accompanying this package, the elrm package vignette may be downloaded from the elrm package homepage on the SFU Statistical Genetics Working Group website (http://stat-db.stat.sfu.ca:8080/statgen/research/elrm). The vignette is also distributed with the code.

References

Forster J.J., McDonald J.W. and Smith P.W.F. Markov Chain Monte Carlo Exact Inference for Binomial and Multinomial Logistic Regression Models. Statistics and Computing, 13:169-177, 2003.

Geyer C.J. Practical Markov chain Monte Carlo. Statistical Science, 7:473-511, 1992.

Zamar David. Monte Carlo Markov Chain Exact Inference for Binomial Regression Models. Master's thesis, Statistics and Actuarial Sciences, Simon Fraser University, 2006.

Zamar D, McNeney B and Graham J. elrm: Software Implementing Exact-like Inference for Logistic Regression Models. Journal of Statistical Software 2007, 21(3).

See Also

update.elrm, summary.elrm, plot.elrm.

Examples

# Drug dataset example with sex and treatment as the variables of interest
data(drugDat);
drug.elrm = elrm(formula=recovered/n~sex+treatment, interest=~sex+treatment, 
	r=4,iter=50000, burnIn=1000, dataset=drugDat);

## Not run: 
# # crash dataset example where the terms of interest are age and 
# # the interaction of age and velocity.
# data(crashDat);
# # The following call produces the error message shown below.
# crash.elrm = elrm(formula=y/n~vel+age+acl+age:vel, interest=~age+age:vel, r=4,
# 	 iter=5000, dataset=crash, burnIn=100);
# 
# # Error in getDesignMatrix(formula, interest, dataset = dataset) : 
# # the 'term.labels' attribute of 'terms.formula(interest)' must match those 
# # found in 'terms.formula(formula)'
# 
# # The error occurs, because the variables within the interaction term in a formula 
# # are re-ordered by the ordering in which the variables occur. Thus, the interaction 
# # between age and velocity is labeled as vel:age in the 'formula' model and as age:vel 
# # in the 'interest' model.
# 
# attr(terms.formula(y/n~vel+age+acl+age:vel),"term.labels");
# # [1] "vel"     "age"     "acl"     "vel:age"
# 
# attr(terms.formula(~age+age:vel),"term.labels");
# # [1] "age"     "age:vel"
# 
# # To get around this problem, place age before vel in the 'formula' model.
# crash.elrm = elrm(formula=y/n~age+vel+acl+age:vel, interest=~age+age:vel, r=4, 
# 	iter=5000, dataset=crash, burnIn=100);
# ## End(Not run)

## Not run: 
# # Urinary tract dataset example with dia as the variable of interest
# data(utiDat); 
# uti.elrm = elrm(uti/n~age+current+dia+oc+pastyr+vi+vic+vicl+vis, interest=~dia,r=4, 
# 	iter=30000,burnIn=1000, dataset=utiDat);## End(Not run)

## Not run: 
# # Titanic dataset example where the variable of interest, class, is coded as a factor
# data(titanDat); 
# titanic.elrm = elrm(surv/n~as.factor(class)+age+sex, interest=~as.factor(class), r=4, 
# 	iter=50000, burnIn=1000, dataset=titanDat);## End(Not run)