Last chance! 50% off unlimited learning
Sale ends in
Estimate threshold generalized linear models, Cox proportional hazards models, and linear mixed models. Supports 14 types of two-phase (one threshold) models and 1 type of three-phase (two thresholds) model.
chngptm (formula.1, formula.2, family, data, type = c("hinge",
"M01", "M02", "M03", "M04", "upperhinge", "M10",
"M20", "M30", "M40", "M21", "M12", "M21c", "M12c",
"M22", "M22c", "M31", "M13", "M33c", "segmented",
"M11", "segmented2", "M111", "step", "stegmented"),
formula.strat = NULL, weights = NULL, offset = NULL,
REML = TRUE, re.choose.by.loglik = FALSE, est.method =
c("default", "fastgrid2", "fastgrid", "grid",
"smoothapprox"), var.type = c("default", "none",
"robust", "model", "bootstrap", "all"), aux.fit =
NULL, lb.quantile = 0.05, ub.quantile = 0.95,
grid.search.max = Inf, test.inv.ci = TRUE,
boot.test.inv.ci = FALSE, bootstrap.type =
c("nonparametric", "wild", "sieve", "wildsieve",
"awb"), m.out.of.n = 0, subsampling = 0, order.max =
10, ci.bootstrap.size = 1000, alpha = 0.05, save.boot
= TRUE, b.transition = Inf, tol = 1e-04, maxit = 100,
chngpt.init = NULL, search.bound = 10, keep.best.fit =
TRUE, ncpus = 1, verbose = FALSE, ...)
chngptm.xy(x, y, type=c("step","hinge","segmented","segmented2","stegmented"),
...)# S3 method for chngptm
coef (object, ...)
# S3 method for chngptm
residuals (object, ...)
# S3 method for chngptm
vcov (object, var.type=NULL, ...)
# S3 method for chngptm
print (x, ...)
# S3 method for chngptm
predict (object, newdata = NULL,
type = c("link", "response", "terms"), ...)
# S3 method for chngptm
plot (x, which = NULL, xlim = NULL, ylim = NULL, lwd = 2,
lcol = "red", lty = 1, add = FALSE, add.points = TRUE,
add.ci = TRUE, breaks = 20, mark.chngpt = TRUE, xlab =
NULL, ylab = NULL, plot.individual.line = FALSE, main
= "", y.adj = NULL, auto.adj.y = FALSE, transform =
NULL, ...)
# S3 method for chngptm
summary (object, var.type = NULL, expo = FALSE,
show.slope.post.threshold = FALSE, verbose = FALSE,
boot.type = "perc", ...)
# S3 method for chngptm
logLik (object, ...)
# S3 method for chngptm
AIC (object, ...)
lincomb(object, comb, alpha = 0.05, boot.type = "perc")
A an object of type chngptm with the following components
Boolean
vector. Estimated coefficients. The last element, named ".chngpt", is the estimated change point
htest. Max score test results
integer. Number of iterations
The part of formula that is free of terms involving thresholded variables
The part of formula that is only composed of thresholded variables
stratification formula
string. coxph or any valid argument that can be passed to glm. But variance estimate is only available for binomial and gaussian (only model-based for latter)
data frame.
type
transform
Numeric. Controls whether threshold model or smooth transition model. Default to Inf, which correponds to threshold model
default: estimation algorithm will be chosen optimally; fastgrid2: a super fast grid search algorithm, limited to linear regression; grid: plain grid search, works for almost all models; smoothapprox: approximates the likelihood function using a smooth function, only works for some models. fastgrid = fastgrid2, kept for backward compatibility
string. Different methods for estimating covariance matrix and constructing confidence intervals
a model fit object that is needed for model-robust estimation of covariance matrix
The maximum number of grid points used in grid search. When doing fast grid search, grid.search.max is set to Inf internally because it does not take more time to examine all potential thresholds.
Boolean, whether or not to find test-inversion confidence interval for threshold
integer, number of bootstrap
double, norminal type I error rate
Boolean, whether to save bootstrap samples
lower bound of the search range for change point estimate
upper bound of the search range for change point estimate
Numeric. Stopping criterion on the coefficient estimate.
integer. Maximum number of iterations in the outer loop of optimization.
numeric. Initial value for the change point.
passed to glm
Boolean.
Boolean.
Boolean.
Boolean.
integer.
Number of cores to use if the OS is not Windows.
Boolean.
outcome
boolean
chngptm fit object.
newdata
chngptm fit object.
arguments passed to glm or coxph
sample size for m-out-of-n bootstrap, default 0 for not doing this type of bootstrap
sample size for subsampling bootstrap, default 0 for not doing this type of bootstrap
whether to get test inversion CI under bootstrap
bounds for search for sloping parameters
an integer
y.adj
auto.adj.y
xlim
ylim
lwd
line col
mark.chngpt
xlab
ylab
offset
lty
lty
nonparametric: the default, classical Efron bootstrap, works for homoscedastic and heteroscedastic indepdendent errors; sieve: works for homoscedastic autocorrelated errors; wild: works for heteroscedastic independent errors; wildsieve: works for heteroscedastic autocorrelated errors; awb: autoregressive wild bootstrap, also works for heteroscedastic autocorrelated errors, but performance may not be as good as wildsieve
order of autocorrelation for autocorrelated errors in sieve and wildsieve bootstrap
a vector of combination coefficients that will be used to form an inner product with the estimated slope
If family is binomial and expo is TRUE, coefficients summary will be shown on the scale of odds ratio instead of slopes
mixed model fitting - should the estimates be chosen to optimize the REML criterion for a fixed threshold
mixed model fitting - should the estimates be chosen to optimize likelihood (REML nor not) or goodness of fit
boolean
character string
Without lb.quantile and ub.quantile, finite sample performance of estimator drops considerably!
When est.method is smoothapprox, Newton-Raphson is done with initial values chosen by change point hypothesis testing. The testing procedure may be less subjective to finite sample volatility.
If var.method is bootstrap, summary of fitted model contains p values for each estimated slope. These p values are approximate p-values, obtained assuming that the bootstrap distributions are normal.
When var.method is bootstrap and the OS is not Windows, the boot package we use under the hood takes advantage of ncpus cores through parallel::mclapply.
lincomb can be used to get the estimate and CI for a linear combination of slopes.
Son, H, Fong, Y. (2020) Fast Grid Search and Bootstrap-based Inference for Continuous Two-phase Polynomial Regression Models, Environmetrics, in press.
Elder, A., Fong, Y. (2020) Estimation and Inference for Upper Hinge Regression Models, Environmental and Ecological Statistics, 26(4):287-302.
Fong, Y. (2019) Fast bootstrap confidence intervals for continuous threshold linear regression, Journal of Computational and Graphical Statistics, 28(2):466-470.
Fong, Y., Huang, Y., Gilbert, P., Permar S. (2017) chngpt: threshold regression model estimation and inference, BMC Bioinformatics, 18(1):454.
Fong, Y., Di, C., Huang, Y., Gilbert, P. (2017) Model-robust inference for continuous threshold regression models, Biometrics, 73(2):452-462.
Pastor-Barriuso, R. and Guallar, E. and Coresh, J. (2003) Transition models for change-point estimation in logistic regression. Statistics in Medicine. 22:13141
# also see the vignette for examples
# threshold linear regression
# for actual use, set ci.bootstrap.size to default or higher
par(mfrow=c(2,2))
types=c("hinge", "segmented", "M02", "M03")
for (type in types) {
fit=chngptm(formula.1=logratio~1, formula.2=~range, lidar, type=type, family="gaussian",
var.type="bootstrap", ci.bootstrap.size=100)
print(summary(fit))
for (i in 1:3) plot(fit, which=i)
out=predict(fit)
plot(lidar$range, out, main=type)
}
# with weights
dat.1=sim.chngpt("thresholded", "segmented", n=200, seed=1, beta=1, alpha=-1, x.distr="norm", e.=4,
family="gaussian")
fit.1.a=chngptm(formula.1=y~z, formula.2=~x, family="gaussian", dat.1, type="segmented",
est.method="fastgrid", var.type="bootstrap", weights=ifelse(dat.1$x<3.5,100,1)
, ci.bootstrap.size=10)
summary(fit.1.a)
plot(fit.1.a)
# fit.1.a$vcov$boot.samples
if (FALSE) {
# likelihood test, combination of slopes
dat=sim.chngpt("thresholded", "segmented", n=200, seed=1, beta=1, alpha=-1, x.distr="norm", e.=4,
family="gaussian")
fit=chngptm(y~z, ~x, family="gaussian", dat, type="segmented", ci.bootstrap.size=100)
fit.0=lm(y~1,dat)
# likelihood ratio test using lmtest::lrtest
library(lmtest)
lrtest(fit, fit.0)
# estimate the slope after threshold using lincomb function in the chngpt package
lincomb(fit, c(0,0,1,1))
}
# threshold logistic regression
dat.2=sim.chngpt("thresholded", "step", n=200, seed=1, beta=1, alpha=-1, x.distr="norm", e.=4,
family="binomial")
fit.2=chngptm(formula.1=y~z, formula.2=~x, family="binomial", dat.2, type="step", est.method="grid")
summary(fit.2)
# no variance estimates available for discontinuous threshold models such as step
# vcov(fit.2$best.fit) gives the variance estimates for the best model conditional on threshold est
# also supports cbind() formula on left hand side
set.seed(1)
dat.2$success=rbinom(nrow(dat.2), 10, 1/(1 + exp(-dat.2$eta)))
dat.2$failure=10-dat.2$success
fit.2a=chngptm(formula.1=cbind(success,failure)~z, formula.2=~x, family="binomial", dat.2,
type="step")
# Poisson example
counts <- c(18,17,15,20,10,20,25,13,12,33,35)
x <- 1:length(counts)
print(d.AD <- data.frame(x, counts))
fit.4=chngptm(formula.1=counts ~ 1, formula.2=~x, data=d.AD, family="poisson",
type="segmented", var.type="bootstrap", verbose=1, ci.bootstrap.size=1)
summary(fit.4)
fit.4a=chngptm(formula.1=counts ~ 1, formula.2=~x, data=d.AD, family="quasipoisson",
type="segmented", var.type="bootstrap", verbose=1, ci.bootstrap.size=1)
if (FALSE) {
# Not run because otherwise the examples take >5s and that is a problem for R CMD check
# coxph example
library(survival)
fit=chngptm(formula.1=Surv(time, status) ~ ph.ecog, formula.2=~age, data=lung, family="coxph",
type="segmented", var.type="bootstrap", ci.bootstrap.size=10)
summary(fit)
# one interaction term (mtcars is part of R default installation)
# est.method will be grid as fastgrid not available for models with interaction terms yet
fit=chngptm(formula.1=mpg ~ hp, formula.2=~hp*drat, mtcars, type="segmented",
family="gaussian", var.type="bootstrap", ci.bootstrap.size=10)
summary(fit)
# interaction, upperhinge model, bootstrap
fit=chngptm(formula.1=mpg ~ hp, formula.2=~hp*drat, mtcars, type="M10",
family="gaussian", var.type="bootstrap", ci.bootstrap.size=10)
summary(fit)
# more than one interaction term
# subsampling bootstrap confidence interval for step model
fit=chngptm(formula.1=mpg~hp+wt, formula.2=~hp*drat+wt*drat, mtcars, type="step",
family="gaussian", var.type="bootstrap", ci.bootstrap.size=10)
summary(fit)
# step model, subsampling bootstrap confidence intervals
fit=chngptm(formula.1=mpg~hp, formula.2=~drat, mtcars, type="step",
family="gaussian", var.type="bootstrap", ci.bootstrap.size=10, verbose=TRUE)
summary(fit)
# higher order threshold models
dat=sim.chngpt(mean.model="thresholded", threshold.type="M22", n=500, seed=1,
beta=c(32,2,10, 10), x.distr="norm", e.=6, b.transition=Inf, family="gaussian",
alpha=0, sd=0, coef.z=0)
fit.0=chngptm(formula.1=y~z, formula.2=~x, dat, type="M22", family="gaussian",
est.method="fastgrid2"); plot(fit.0)
dat=sim.chngpt(mean.model="thresholded", threshold.type="M22c", n=500, seed=1,
beta=c(32,2,32, 10), x.distr="norm", e.=6, b.transition=Inf, family="gaussian",
alpha=0, sd=0, coef.z=0)
fit.0=chngptm(formula.1=y~z, formula.2=~x, dat, type="M22c", family="gaussian",
est.method="fastgrid2"); plot(fit.0)
# examples of aux.fit
fit.0=glm(yy~zz+ns(xx,df=3), data, family="binomial")
fit = chngptm (formula.1=yy~zz, formula.2=~xx, family="binomial", data, type="hinge",
est.method="smoothapprox", var.type="all", verbose=verbose, aux.fit=fit.0,
lb.quantile=0.1, ub.quantile=0.9, tol=1e-4, maxit=1e3)
}
# example of random intercept
dat=sim.twophase.ran.inte(threshold.type="segmented", n=50, seed=1)
fit = chngptm (formula.1=y~z+(1|id), formula.2=~x, family="gaussian", dat,
type="segmented", est.method="grid", var.type="bootstrap", ci.bootstrap.size=1)
plot(fit)
out=predict(fit, re.form=NA)
plot(dat$x, out)
out.1=predict(fit, type="response", re.form=NULL)# includes re
plot(dat$x, out.1, type="p", xlab="x")
Run the code above in your browser using DataLab