Type: | Package |
Title: | Repeated Measures Proportional Odds Logistic Regression |
Version: | 3.4 |
Date: | 2016-02-26 |
Author: | Nick Parsons |
Maintainer: | Nick Parsons <nick.parsons@warwick.ac.uk> |
Description: | Fits linear models to repeated ordinal scores using GEE methodology. |
License: | GPL-3 |
Imports: | Rcpp (≥ 0.11.3), Matrix, methods |
LinkingTo: | Rcpp, RcppArmadillo |
NeedsCompilation: | yes |
Packaged: | 2016-02-26 08:36:36 UTC; mlsgbf |
Repository: | CRAN |
Date/Publication: | 2016-02-27 10:11:34 |
Repeated Measures Proportional Odds Logistic Regression using GEE
Description
The package allows regression models to be fitted to repeated ordinal scores, for the proportional odds model, using a generalized estimating equation (GEE) methodology. The algorithm estimates the correlation parameter by minimizing the generalized variance of the regression parameters at each step of the fitting algorithm. Parameter estimation is available for the uniform and first-order autoregressive correlation models, for data potentially recorded at irregularly spaced time intervals. A test for the proportional odds assumption is available, as is an option to fit a polynomial model to the the cut-point parameters.
Usage
repolr(formula, subjects, data, times, categories, corr.mod = "independence",
alpha = 0.5, po.test = FALSE, fixed = FALSE, poly = NULL,
space = NULL, diffmeth = "analytic", fit.opt = rep(NA, 5))
Arguments
formula |
a formula, as for other regression models. |
subjects |
a character string specifying the name of the subject variable. |
data |
a dataframe in which to interpret the variables occurring in the formula. |
times |
a vector of times which occur within subject clusters; e.g. for four evenly spaced times c(1, 2, 3, 4). |
categories |
a numeric variable indicating the number of ordinal score categories. |
corr.mod |
a character string specifying the correlation structure. The following are permitted: “ |
alpha |
an initial value for the correlation parameter. |
po.test |
a logical variable; if true a score test for proportional odds is reported. |
fixed |
a logical variable; if true the correlation is fixed at the initial value ( |
poly |
a numeric variable indicating the order of the polynomial contrasts used for the cut-point model. |
space |
a vector indicating the category spacing when fitting the polynomial model; can generally be set to 1: |
diffmeth |
a character string specifying the method used for estimation of |
fit.opt |
a vector of options to control the behaviour of the fitting algorithm. |
Details
The repolr
function fits models for repeated ordinal scores using GEE methodology.
The user is required to specify, as a minimum: (i) a data set name (data
), (ii) a model formula (formula
), (iii) a cluster identification variable (subjects
), (iv) a time variable (time
) and (v) the number of categories used for the response variable (categories
).
The data may contain records with missing data for either the response variable or the explanatory variables. The response variable must have at least three ordered categories (K greater than or equal to 3) and, as K-1 cut-point parameters are estimated, an additional intercept term can not be explicitly included in the model formula. A subject variable, which takes integer values from 1 to N with no missing values allowed, indicates the data clusters (patients or experimental units) and a time variable indicates the within cluster ordering; times must be ordered integers starting from one and spaced to indicate the relative distance between successive times. For instance, four observations at equally spaced times would be entered as 1, 2, 3 and 4, whereas if the first two observations were separated by half the time interval of the other observations then coding would be 1, 2, 4 and 6. The data must be pre-sorted by time clusters within each subject, and complete, i.e. where data is missing for a particular time it must be indicated as such. the datasets provided with this package provide exemplars of the required data formatting; e.g. HHSpain
and QoL
.
The available options for the correlation model (corstr
) are AR1, uniform, fixed and independence, with default setting independence.
Additionally there are a number of other algorithm related options.
The algorithm is generally robust to the initial value for alpha
(default setting = 0.5), where estimation is rerquired, however a starting value for alpha
can be set. If required the correlation parameter, set via alpha
, can be fixed throughout model fitting, and not updated, by setting the option fixed
to TRUE
.
The partial derivatives of the log of the determinant of the robust variance matrix (generalized variance), with respect to alpha
, can either be determined analytically (“analytic
” setting for diffmeth
) or numerically by finite differencing (“numeric
” setting for diffmeth
). The latter method is often quicker for complex regression models, or if K is large.
Function poly
, specifies the order of fitted orthogonal polynomial contrasts for the cut-point parameters; the default setting for repolr is to fit the complete set of cut-point parameters. Fitting polynomial contrasts can be particularly useful for long ordinal scores, where K is large, or where a particular form for the cut-point parameters is preferred; e.g. increasing uniformly at the extremes of the score scale. The order of the polynomial must be an integer less than K-1. The function requires one additional argument, space
, that indicates the spacing between categories. This is normally set to 1:categories
when all categories are observed, but can be modified to provide more realsitic models if for instance some categories are not observed.
The proportional odds assumption can be tested, using a score test, by setting the option po.test
to TRUE.
The fit.opt
argument provides control of the fitting algorithm; the defaults are
c(cmaxit = 10, omaxit = 5, ctol = 0.001, otol = 0.00001, h = 0.01)
.
These control the maximum number of iterations for updating estimates of alpha
, the maximum number of iterations for updating the regression parameters within each of the updating steps for alpha
, the convergence tolerances for estimation of alpha
and the regression parameters, and the interval h
for finite differencing, if the “numeric
” option is selected.
Model fitting is implemented via a suite of of functions developed using RcppArmadillo (Rcpp) that construct correlations between derived binary variables at each time-point (smat
) and between time-points (cmat
). Complete (sparse) model covariance matrices are constructed using hgmat
and alphpow
, with model parameter estimation implemented in ordgee
, using current estimates of the inverse of the correlation matrix from icormat
. Function upalpha
provides updates of correlation parameter estimates, and potest
implements the test of proportional odds. These functions are not documented in detail here as they are primarily for internal use within repolr. There use outside of this setting is not recommended.
Value
The function summary.repolr
is used to obtain and print a summary of the fitted model.
The fitted model is an object of class “repolr
” and has the following values:
poly.mod |
polynomial model for cut points: a |
y |
the response variable. |
linear.predictors |
a vector of linear predictors. |
fitted.values |
a vector of the fitted values. |
coefficients |
a named vector of regression coefficients. |
robust.var |
the robust (sandwich) variance matrix. |
naive.var |
the naive variance matrix. |
alpha |
an estimate of the correlation parameter. |
convergence |
a logical variable to reported if convergence was achieved. |
iter |
the number of iterations. |
grad1 |
first derivative of generalized variance at convergence. |
grad2 |
second derivative of generalized variance at convergence. |
crit |
convergence criterion. |
po.test |
results of po test: a |
References
Parsons NR, Costa ML, Achten J, Stallard N. Repeated measures proportional odds logistic regression analysis of ordinal score data in the statistical software package R. Computational Statistics and Data Analysis 2009; 53:632-641.
Parsons NR, Edmondson RN, Gilmour SG. A generalized estimating equation method for fitting autocorrelated ordinal score data with an application in horticultural research. Journal of the Royal Statistical Society C 2006; 55:507-524.
Stiger TR, Barnhart HX, Williamson JM. Testing proportionality in the proportional odds model fitted with GEE. Statistics in Medicine 1999; 18:1419-1433.
Parsons NR. Proportional-odds models for repeated composite and long ordinal outcome scales. Statistics in Medicine 2013; 32:3181-3191.
See Also
Examples
data(HHSpain)
mod.0 <- repolr(HHSpain~Sex*Time, data=HHSpain, categories=4, subjects="Patient",
times=c(1,2,5), corr.mod="uniform", alpha=0.5)
summary(mod.0)
summary(update(mod.0, diffmeth = "numeric"))
summary(update(mod.0, fixed = TRUE, alpha =0.5))
Harris Hip Pain Scores
Description
Harris hip pain scores, measured on a four point ordinal score scale, post-surgery for 58 patients at three occasions.
Usage
data(HHSpain)
Format
A data frame with 174 observations on the following four variables.
Patient
a patient identifier variable.
Sex
a factor with levels
F
orM
Time
recorded at baseline (1), two years (2) and five years (5) post-surgery.
HHSpain
a pain score coded as none (1), slight (2), mild (3) and moderate or marked (4).
Examples
data(HHSpain)
head(HHSpain)
Quasilikelihood Information Criterion
Description
The quasilikelihood information criterion (QIC) developed by Pan (2001) is a modification of the Akaike information criterion (AIC) for models fitted by GEE. QIC
is used for choosing the best correaltion structure and QICu
is used for choosing the best subset of covariates. The quasilikelihood (QLike
) is also reported for completeness. When choosing between two or more models, with different subset of covariates, the one with the smallest QICu
measure is preferred and similarly, when choosing between competing correlation structures, with the same subset of covariates in both, the model with the smallest QIC
measure is preferred.
Usage
QIC(object, digits = 3)
Arguments
object |
is a fitted model using |
digits |
the number of decimal places to display in reported summaries. |
Value
QLike |
model quasilikelihood. |
QIC |
model |
QICu |
model |
References
Pan W. Akaikes information criterion in generalized estimating equations. Biometrics 2001; 57:120-125.
Examples
data(HHSpain)
mod.0 <- repolr(HHSpain~Time, data=HHSpain, categories=4, subjects="Patient",
times=c(1,2,5), corr.mod="independence", alpha=0.5)
QIC(mod.0)
QIC(update(mod.0, formula = HHSpain~Time + Sex))$QICu
QIC(update(mod.0, formula = HHSpain~Time * Sex))$QICu
Quality of Life Scores
Description
Quality of life scores after hip replacement.
Usage
data(QoL)
Format
A data frame with 336 observations on the following four variables.
QoL
a numeric vector.
Patient
a numeric vector.
Time
a numeric vector.
Treat
a factor with levels
A
B
.
Examples
data(QoL)
head(QoL)
Achilles Tendon Rupture
Description
EuroQol activity scores, measured on a three point ordinal score scale, post-surgery for 48 patients at three occasions.
Usage
data(achilles)
Format
A data frame with 144 observations on the following four variables.
Patient
a patient identifier variable.
Treat
post-surgery treatments are either immediate mobilisation in a carbon-fibre orthosis with three 1.5cm heel raises (1) or traditional plaster cast immobilisation (2).
Time
recorded at baseline (1), six months (2) and one year (3) post-surgery.
Activity
ability to undertake usual activities post-surgery; this was scored by each patient as either no problem (1), some problem (2) or an inability (3) to perform usual activity (e.g. work, leisure, housework etc).
Source
Costa, M.L., MacMillan, K., Halliday, D., Chester, R., Shepstone, L., Robinson, A.H.N., Donell, S.T. (2006). Randomised controlled trials of immediate weight-bearing mobilisation for rupture of the tendon Achillis. Journal of Bone and Joint Surgery (British) 88-B, 69-77.
Examples
data(achilles)
head(achilles)
Begonia Pot Plant Quality Scores
Description
Begonia pot plant quality scores (for 2 varieties and 3 transport chains), during 5 weeks in simulated shelf-life conditions (temperature, light and irrigation). Quality scores were originally made on an ordinal scale from 1 to 10 (highest quality). However, only categories 2 to 9 were used, so these were re-coded to scale from 1 to 8. In addition to overall quality scores, a range of plant physiological characteristics were also observed.
Usage
data(begonia)
Format
A data frame with 720 observations on the following variables.
Pot
pot plant identifier; 1 to 144.
Plot
location in growing compartment; 1 to 48.
Week
week number in simulated shelf-life; 1 to 5.
Temp
temperature in compartment;
16C
or21C
.Light
light level in compartment;
High
orLow
.Chain
transport chain;
Cold
,Comm
(commercial) orOptm
(optimum).Irrig
irrigation in compartment;
Cont
(control) orFluct
(fluctuating).Variety
variety;
Balli
orBatik
.Qual
quality score; 1 to 8.
FDrop
count of dropped flowers.
CBDrop
count of dropped coloured buds.
GBDrop
count of dropped green buds.
FDam
count of damaged flowers.
FSing
count of single flowers.
FDoub
count of double flowers.
LYell
count of yellow leaves.
Source
Parsons NR, Edmondson RN, Gilmour SG. A generalized estimating equation method for fitting autocorrelated ordinal score data with an application in horticultural research. Journal of the Royal Statistical Society C 2006; 55:507-524.
Examples
data(begonia)
head(begonia)
Confidence Intervals for repolr Model Parameters
Description
Computes confidence intervals for one or more parameters in a fitted repolr
model object.
Usage
## S3 method for class 'repolr'
confint(object, parm, level = 0.95, robust.var = TRUE, ...)
Arguments
object |
is a model fitted using |
parm |
a specification of which parameters are to be used, either a vector of numbers or a vector of names. If missing, all parameters are considered. |
level |
the confidence level required. |
robust.var |
logical; if |
... |
further arguments passed to or from other methods. |
Details
The method assumes normality and uses as default the estimated robust variance matrix.
Value
A matrix (or vector) with columns giving lower and upper confidence limits for each parameter.
Examples
data(HHSpain)
mod.0 <- repolr(HHSpain~Sex*Time, data=HHSpain, categories=4, subjects="Patient",
times=c(1,2,5), corr.mod="uniform", alpha=0.5)
confint(mod.0, robust.var = FALSE)
Mobility after Hip Fracture Fixation Surgery
Description
Patient mobility score after hip fracture fixation surgery, using two distinct procedures.
Usage
data(mobility)
Format
A data frame with 600 observations on the following variables.
subject
patient identifier.
time
assessemnt occasions; month 1, 2, 3 and 4.
treat
intervention group;
A
orB
.age
patient age in years.
gender
gender;
F
orM
.mobility
mobility score; 1, 2 and 3.
Examples
data(mobility)
head(mobility)
Expand Ordinal Data for Model Fitting with repolr
Description
Expands ordinal score data into an appropriate form, K-1 new binary variables for an ordinal score with K categories, for model fitting using repolr
.
Usage
ord.expand(space, formula, times, poly, data, subjects, categories)
Arguments
space |
a vector indicating the category spacing when fitting the polynomial model. |
formula |
a formula, as for other regression models. |
times |
a vector of times which occur within subject clusters. |
poly |
a numeric variable indicating the order of the polynomial contrasts used for the cut-point model. |
data |
a dataframe in which to interpret the variables occurring in the formula. |
subjects |
a character string specifying the name of the subject variable. |
categories |
a numeric variable indicating the number of ordinal score categories. |
Details
For internal use with repolr
.
Value
An expanded dataframe.
Estimates Cut-point Parameters for Fitted repolr Model
Description
After fitting a model using repolr
, function polycuts
gives estimates and standard errors for the K-1 cut-point parameters, based on the polynomial model from the fit of repolr
. Polynomial cut-point parameter estimates from the orginal model are also shown.
Usage
polycuts(object, digits = 3, robust.var = TRUE)
Arguments
object |
is a model fitted using |
digits |
the number of decimal places to display in reported summaries. |
robust.var |
a logical variable: if |
Value
coef |
polynomial parameter estimates from |
poly |
a vector of K-1 cut-point parameters. |
order |
the order of the polynomial. |
Examples
data(HHSpain)
mod.0 <- repolr(HHSpain~Sex*Time, data=HHSpain, categories=4, subjects="Patient",
times=c(1,2,5), corr.mod="uniform", alpha=0.5)
summary(mod.0)
mod.1 <- update(mod.0, poly=1)
summary(mod.1)
polycuts(mod.1)
mod.2 <- update(mod.0, poly=2)
summary(mod.2)
polycuts(mod.2)
Predict Method for Fitted repolr Model
Description
Calculates predictions and standard errors of predictions for a fitted repolr
model object.
Usage
## S3 method for class 'repolr'
predict(object, newdata = NULL, se.fit = FALSE,
robust.var = TRUE, type = c("link", "response", "terms"), ...)
Arguments
object |
is a model fitted using |
newdata |
optionally, a data frame in which to find variables with which to predict; if missing the model fitted values are reported. |
se.fit |
Logical indicating if standard errors are required. |
robust.var |
logical; if |
type |
is the type of prediction required. The default “ |
... |
further arguments passed to or from other methods. |
Details
If newdata is missing predictions are based on the data used to fit the repolr
model. If newdata are supplied then the format of these data must conform to the same format required for model fitting using repolr
. See repolr
for details.
Value
fit |
Predictions. |
se.fit |
Estimated standard errors. |
Examples
data(HHSpain)
mod.0 <- repolr(HHSpain~Sex*Time, data=HHSpain, categories=4, subjects="Patient",
times=c(1,2,5), corr.mod="uniform", alpha=0.5)
predict(mod.0, newdata = data.frame(Patient = rep(100, 3), Time = c(1, 2, 5),
Sex = factor(rep(1, 3), levels=1:2, labels=c("F", "M"))),
type="link", se.fit = TRUE)
Summary of Fitted repolr Model
Description
Function to summarise the fit of a repolr model.
Usage
## S3 method for class 'repolr'
summary(object, digits, robust.var, ...)
Arguments
object |
fitted model. |
digits |
integer for number formatting. |
robust.var |
logical; if |
... |
further arguments passed to or from other methods. |
Details
Default is to use robust variance estimates. However, if robust.var
is set to FALSE
, naive variance estimates are used.
Calculates the Variance-Covariance Matrix for Fitted repolr Model
Description
Returns the variance-covariance matrix of the main parameters of a fitted repolr
model object.
Usage
## S3 method for class 'repolr'
vcov(object, robust.var = TRUE, ...)
Arguments
object |
is a model fitted using |
robust.var |
logical; if |
... |
further arguments passed to or from other methods. |
Details
Default is to output the estimated robust variance matrix. However, if robust.var
is set to FALSE
, the naive variance matrix is reported.
Value
A matrix of the estimated covariances between the parameter estimates.
Examples
data(HHSpain)
mod.0 <- repolr(HHSpain~Sex*Time, data=HHSpain, categories=4, subjects="Patient",
times=c(1,2,5), corr.mod="uniform", alpha=0.5)
vcov(mod.0, robust.var = FALSE)
Working Correlation Matrix for Fitted repolr Model
Description
The working correlation matrix for the selected model; “ar1
”, “uniform
” or “independence
”.
Usage
work.corr(object, digits = 3)
Arguments
object |
is a model fitted using |
digits |
integer for number formatting. |
Value
A T(K-1) correlation matrix.
Examples
data(HHSpain)
mod.0 <- repolr(HHSpain~Sex*Time, data=HHSpain, categories=4, subjects="Patient",
times=c(1,2,5), corr.mod="uniform", alpha=0.5)
work.corr(mod.0)