chngptm: Threshold Models Estimation

In chngpt: Estimation and Hypothesis Testing for Threshold Regression

Threshold Models Estimation

Description

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.

Usage

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 class 'chngptm'
 coef(object, ...)
## S3 method for class 'chngptm'
 residuals(object, ...)
## S3 method for class 'chngptm'
 vcov(object, var.type=NULL, ...)
## S3 method for class 'chngptm'
 print(x, ...)
## S3 method for class 'chngptm'
 predict(object, newdata = NULL, 
 type = c("link", "response", "terms"), ...)
## S3 method for class '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 class 'chngptm'
 summary(object, var.type = NULL, expo = FALSE,
 show.slope.post.threshold = FALSE, verbose = FALSE,
 boot.type = "perc", ...)
## S3 method for class 'chngptm'
 logLik(object, ...)
## S3 method for class 'chngptm'
 AIC(object, ...)
 

lincomb(object, comb, alpha = 0.05, boot.type = "perc")

Arguments

formula.1	The part of formula that is free of terms involving thresholded variables
formula.2	The part of formula that is only composed of thresholded variables
formula.strat	stratification formula
family	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	data frame.
type	type
transform	transform
b.transition	Numeric. Controls whether threshold model or smooth transition model. Default to Inf, which correponds to threshold model
est.method	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
var.type	string. Different methods for estimating covariance matrix and constructing confidence intervals
aux.fit	a model fit object that is needed for model-robust estimation of covariance matrix
grid.search.max	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.
test.inv.ci	Boolean, whether or not to find test-inversion confidence interval for threshold
ci.bootstrap.size	integer, number of bootstrap
alpha	double, norminal type I error rate
save.boot	Boolean, whether to save bootstrap samples
lb.quantile	lower bound of the search range for change point estimate
ub.quantile	upper bound of the search range for change point estimate
tol	Numeric. Stopping criterion on the coefficient estimate.
maxit	integer. Maximum number of iterations in the outer loop of optimization.
chngpt.init	numeric. Initial value for the change point.
weights	passed to glm
verbose	Boolean.
add.points	Boolean.
add.ci	Boolean.
add	Boolean.
breaks	integer.
ncpus	Number of cores to use if the OS is not Windows.
keep.best.fit	Boolean.
y	outcome
show.slope.post.threshold	boolean
x	chngptm fit object.
newdata	newdata
object	chngptm fit object.
...	arguments passed to glm or coxph
m.out.of.n	sample size for m-out-of-n bootstrap, default 0 for not doing this type of bootstrap
subsampling	sample size for subsampling bootstrap, default 0 for not doing this type of bootstrap
boot.test.inv.ci	whether to get test inversion CI under bootstrap
search.bound	bounds for search for sloping parameters
which	an integer
y.adj	y.adj
auto.adj.y	auto.adj.y
xlim	xlim
ylim	ylim
lwd	lwd
lcol	line col
mark.chngpt	mark.chngpt
xlab	xlab
ylab	ylab
offset	offset
lty	lty
boot.type	lty
bootstrap.type	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.max	order of autocorrelation for autocorrelated errors in sieve and wildsieve bootstrap
comb	a vector of combination coefficients that will be used to form an inner product with the estimated slope
expo	If family is binomial and expo is TRUE, coefficients summary will be shown on the scale of odds ratio instead of slopes
REML	mixed model fitting - should the estimates be chosen to optimize the REML criterion for a fixed threshold
re.choose.by.loglik	mixed model fitting - should the estimates be chosen to optimize likelihood (REML nor not) or goodness of fit
plot.individual.line	boolean
main	character string

Details

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.

Value

A an object of type chngptm with the following components

converged	Boolean
coefficients	vector. Estimated coefficients. The last element, named ".chngpt", is the estimated change point
test	htest. Max score test results
iter	integer. Number of iterations

References

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

Examples

# 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

## Not run: 
# 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))

## End(Not run)


# 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)



## Not run: 
# 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)




## End(Not run)

# 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")

chngpt documentation built on Feb. 16, 2023, 10:38 p.m.