Package 'MultiKink' reference manual

Title:	Estimation and Inference for Multi-Kink Quantile Regression
Description:	Estimation and inference for multiple kink quantile regression for longitudinal data and the i.i.d data. A bootstrap restarting iterative segmented quantile algorithm is proposed to estimate the multiple kink quantile regression model conditional on a given number of change points. The number of kinks is also allowed to be unknown. In such case, the backward elimination algorithm and the bootstrap restarting iterative segmented quantile algorithm are combined to select the number of change points based on a quantile BIC. For longitudinal data, we also develop the GEE estimator to incorporate the within-subject correlations. A score-type based test statistic is also developed for testing the existence of kink effect. The package is based on the paper, ``Wei Zhong, Chuang Wan and Wenyang Zhang (2022). Estimation and inference for multikink quantile regression, JBES'' and ``Chuang Wan, Wei Zhong, Wenyang Zhang and Changliang Zou (2022). Multi-kink quantile regression for longitudinal data with application to progesterone data analysis, Biometrics".
Authors:	Chuang Wan [aut, cre], Wei Zhong [aut]
Maintainer:	Chuang Wan <[email protected]>
License:	GPL
Version:	0.2.0
Built:	2024-10-30 03:29:42 UTC
Source:	https://github.com/cran/MultiKink

Auxiliary parameters to control the model fitting

Description

This function defines auxiliary parameters that control the model fitting process.

Usage

fit.control(
  toll = 1e-04,
  h = 1,
  it.max = 50,
  K.max = 6,
  stop.if.error = TRUE,
  dev0 = NULL,
  visual = FALSE,
  visualBoot = FALSE,
  pow = c(1, 1),
  digits = NULL,
  grid = NULL,
  n.boot = 20
)
fit.control(
  toll = 1e-04,
  h = 1,
  it.max = 50,
  K.max = 6,
  stop.if.error = TRUE,
  dev0 = NULL,
  visual = FALSE,
  visualBoot = FALSE,
  pow = c(1, 1),
  digits = NULL,
  grid = NULL,
  n.boot = 20
)

Arguments

`toll`	Positive convergence tolerance.
`h`	Positive factor (from zero to one) modifying the increments in kink parameter updates during the iterative process.
`it.max`	Positive integer for the maximal number of iterations.
`K.max`	Positive integer for the maximal given number of kink points.
`stop.if.error`	Logical indicating if the estimation algorithm should be stopped if some kink point estimators belong to the non-admissible set. Default is FALSE which suggests removing the non-admissible change points automatically.
`dev0`	Initial objective value or deviance. Default is NULL which implies that the initial value is unknown.
`visual`	Logical indicating if the results of the estimation process should be printed at each iteration.
`visualBoot`	Logical indicating if the results of estimation should be printed at each iteration in the bootstrap restarting process.
`pow`	The powers of the pseudo covariates employed by the algorithm.
`digits`	If specified, it means the desired number of decimal points of the kink estimators to be used during the iterative algorithm.
`grid`	It measures how close between the two adjacent change points should be merged, default is NULL.
`n.boot`	Positive integer indicating the times of bootstrap re-sampling in the bootstrap restarting algorithm, default is 20.

Value

A list with the arguments as components to be used by mkqr.fit and mkqr.bea.

Examples

# Example usage
fit.control(K.max=8)

# Example usage
fit.control(K.max=8)

Test the existence of kink effect in the multi-kink quantile regression

Description

This function tests the existence of a kink effect in the multi-kink quantile regression.

Usage

KinkTest(
  y,
  thre.x,
  cont.z,
  id,
  tau = 0.5,
  NB = 200,
  sparsity = "nid",
  bandwidth_type = c("Hall-Sheather", "Bofinger", "Chamberlain")
)
KinkTest(
  y,
  thre.x,
  cont.z,
  id,
  tau = 0.5,
  NB = 200,
  sparsity = "nid",
  bandwidth_type = c("Hall-Sheather", "Bofinger", "Chamberlain")
)

Arguments

`y`	A numeric vector of response.
`thre.x`	A numeric vector of scalar covariate with threshold effect.
`cont.z`	A numeric matrix of design covariates with constant slopes.
`id`	A numeric vector of index used for longitudinal data; can be missing or NULL for iid data.
`tau`	A numeric scalar representing the quantile level (default is 0.5).
`NB`	An integer specifying the resampling times (default is 200).
`sparsity`	The error term type. Specify one from "iid" and "nid" (default is "nid").
`bandwidth_type`	The bandwidth type. Specify one from "Hall-Sheather", "Bofinger", or "Chamberlain" (default is "Hall-Sheather").

Value

A list containing the p-value (pv), the statistic based on the original data (Tn), and the statistics by wild bootstrap (Tn.NB).

Examples

# Example 1: i.i.d data type
library(quantreg)
n = 200
Z1 <- rexp(n,1)
Z2 <- rbinom(n,1,0.5)
Z <- cbind(Z1,Z2)
epsilon <- rnorm(n,0,1)
X <- runif(n,-2,1)
psi <- c(-1,0)
k <- length(psi)
PSI <- matrix(rep(psi,rep(n,k)),ncol=k)
XP <- matrix(rep(X,k),nrow=n)
XR <- cbind(1,X,pmax((XP-PSI),0),Z)
bet <- c(1,-1,0,0,sqrt(3),-sqrt(3))
Y <- XR %*% bet + epsilon
obj <- KinkTest(y=Y,thre.x=X,cont.z=Z,
                bandwidth_type=c("Hall-Sheather"))
obj$pv

## Not run: 
# Example 2: longitudinal data
library(quantreg)
N = 200
T = 5
subject = rep(1:N,each=T)
NT = N*T
Z1 <- rexp(NT,1)
Z2 <- rbinom(NT,1,0.5)
Z <- cbind(Z1,Z2)
epsilon <- rnorm(NT,0,1)
X <- runif(NT,-2,1)
psi <- c(-1,0)
k <- length(psi)
PSI <- matrix(rep(psi,rep(NT,k)),ncol=k)
a <- rnorm(N,0,1)
A <- rep(a,each=T)
XP <- matrix(rep(X,k),nrow=NT)
XR <- cbind(1,X,pmax((XP-PSI),0),Z)
bet <- c(1,-1,0,0,sqrt(3),-sqrt(3))
Y <- XR %*% bet + A + epsilon
obj <- KinkTest(y=Y,thre.x=X,cont.z=Z,id=subject,
                bandwidth_type=c("Hall-Sheather"))
obj$pv

## End(Not run)
# Example 1: i.i.d data type
library(quantreg)
n = 200
Z1 <- rexp(n,1)
Z2 <- rbinom(n,1,0.5)
Z <- cbind(Z1,Z2)
epsilon <- rnorm(n,0,1)
X <- runif(n,-2,1)
psi <- c(-1,0)
k <- length(psi)
PSI <- matrix(rep(psi,rep(n,k)),ncol=k)
XP <- matrix(rep(X,k),nrow=n)
XR <- cbind(1,X,pmax((XP-PSI),0),Z)
bet <- c(1,-1,0,0,sqrt(3),-sqrt(3))
Y <- XR %*% bet + epsilon
obj <- KinkTest(y=Y,thre.x=X,cont.z=Z,
                bandwidth_type=c("Hall-Sheather"))
obj$pv

## Not run: 
# Example 2: longitudinal data
library(quantreg)
N = 200
T = 5
subject = rep(1:N,each=T)
NT = N*T
Z1 <- rexp(NT,1)
Z2 <- rbinom(NT,1,0.5)
Z <- cbind(Z1,Z2)
epsilon <- rnorm(NT,0,1)
X <- runif(NT,-2,1)
psi <- c(-1,0)
k <- length(psi)
PSI <- matrix(rep(psi,rep(NT,k)),ncol=k)
a <- rnorm(N,0,1)
A <- rep(a,each=T)
XP <- matrix(rep(X,k),nrow=NT)
XR <- cbind(1,X,pmax((XP-PSI),0),Z)
bet <- c(1,-1,0,0,sqrt(3),-sqrt(3))
Y <- XR %*% bet + A + epsilon
obj <- KinkTest(y=Y,thre.x=X,cont.z=Z,id=subject,
                bandwidth_type=c("Hall-Sheather"))
obj$pv

## End(Not run)

Fit the multi-kink quantile regression in the absence of the number of change points.

Description

Fit the multi-kink quantile regression in the absence of the number of change points.

Usage

mkqr.bea(
  y,
  thre.x,
  cont.z,
  id,
  tau = 0.5,
  Cn = 1,
  bandwidth_type = "Hall-Sheather",
  control = fit.control(),
  est.type = "WI",
  wi.type = "general",
  wc.type = "cs"
)
mkqr.bea(
  y,
  thre.x,
  cont.z,
  id,
  tau = 0.5,
  Cn = 1,
  bandwidth_type = "Hall-Sheather",
  control = fit.control(),
  est.type = "WI",
  wi.type = "general",
  wc.type = "cs"
)

Arguments

`y`	A numeric vector of response.
`thre.x`	A numeric vector of scalar covariate with threshold effect.
`cont.z`	A numeric matrix of design with constant slopes.
`id`	A numeric vector of index used for longitudinal data; can be missing or NULL for iid data.
`tau`	The quantile level that belongs to (0,1). Default is 0.5.
`Cn`	A positive number corresponding to different types of BIC. Default is 1.
`bandwidth_type`	The bandwidth type. Specify one from "Hall-Sheather", "Bofinger", and "Chamberlain". Default is "Hall-Sheather".
`control`	A list returned by fit.control.
`est.type`	The estimation type for the longitudinal data. Specify one from "WI", "WC", corresponding to the working independence (WI) estimator and working correlation (WC) estimator. Default is "WI".
`wi.type`	If est.type = "WI", then set the error structure of the variance-covariance matrix estimation. Specify one from "Compound", "AR", and "general".
`wc.type`	If est.type = "WC", then set the correlation structure within subject. Specify one from "ar" and "cs". Default is "cs".

Value

A list containing the estimated number of kink points (n.psi), the fitted quantile objective value (rho), estimated regression coefficients with intercept (bet.est), the estimated standard error of the regression coefficients (bet.se), the estimated change points (psi.est), and the estimated standard errors of threshold parameters (psi.se).

Examples

## Not run: 
# Simple examples for iid data type
n <- 500
Z1 <- rexp(n,1)
Z2 <- rbinom(n,1,0.5)
Z <- cbind(Z1,Z2)
epsilon <- rnorm(n,0,1)
X <- runif(n,-2,1)
psi <- c(-1,0)
k <- length(psi)
PSI <- matrix(rep(psi,rep(n,k)),ncol=k)
XP <- matrix(rep(X,k),nrow=n)
XR <- cbind(1,X,pmax((XP-PSI),0),Z)
bet <- c(1,-1,3,-3,sqrt(3),-sqrt(3))
Y <- XR %*% bet + epsilon

# Estimation setting
tau <- 0.5
K.max <- 5
control <- fit.control(K.max = K.max)
Cn <- 1
mkqr.bea(y = Y, thre.x = X, cont.z = Z, tau = tau, Cn = Cn, control = control)

# Simple examples for longitudinal data
N <- 200
T <- 5
subject <- rep(1:N, each = T)
NT <- N * T
Z1 <- rexp(NT, 1)
Z2 <- rbinom(NT, 1, 0.5)
Z <- cbind(Z1, Z2)
epsilon <- rnorm(NT, 0, 1)
X <- runif(NT, -2, 1)
psi <- c(-1, 0)
k <- length(psi)
PSI <- matrix(rep(psi, rep(NT, k)), ncol = k)
a <- rnorm(N, 0, 1)
A <- rep(a, each = T)
XP <- matrix(rep(X, k), nrow = NT)
XR <- cbind(1, X, pmax((XP - PSI), 0), Z)
bet <- c(1, -1, 3, -3, sqrt(3), -sqrt(3))
Y <- XR %*% bet + A + epsilon
tau <- 0.5
k <- 2

# Example 1: the working independence estimator; the error structure is "general"
est.type <- "WI"
wi.type <- "Compound"
tau <- 0.5
K.max <- 5
control <- fit.control(K.max = K.max)
Cn <- 1
mkqr.bea(y = Y, thre.x = X, cont.z = Z, id = subject, tau = tau, Cn = Cn,
         control = control, est.type = est.type, wi.type = wi.type)

# Example 2: the working correlated estimator; the correlation structure is "cs"
est.type <- "WC"
wc.type <- "cs"
tau <- 0.5
K.max <- 5
control <- fit.control(K.max = K.max)
Cn <- 1
mkqr.bea(y = Y, thre.x = X, cont.z = Z, id = subject, tau = tau, Cn = Cn,
         control = control, est.type = est.type, wc.type = wc.type)

## End(Not run)

## Not run: 
# Simple examples for iid data type
n <- 500
Z1 <- rexp(n,1)
Z2 <- rbinom(n,1,0.5)
Z <- cbind(Z1,Z2)
epsilon <- rnorm(n,0,1)
X <- runif(n,-2,1)
psi <- c(-1,0)
k <- length(psi)
PSI <- matrix(rep(psi,rep(n,k)),ncol=k)
XP <- matrix(rep(X,k),nrow=n)
XR <- cbind(1,X,pmax((XP-PSI),0),Z)
bet <- c(1,-1,3,-3,sqrt(3),-sqrt(3))
Y <- XR %*% bet + epsilon

# Estimation setting
tau <- 0.5
K.max <- 5
control <- fit.control(K.max = K.max)
Cn <- 1
mkqr.bea(y = Y, thre.x = X, cont.z = Z, tau = tau, Cn = Cn, control = control)

# Simple examples for longitudinal data
N <- 200
T <- 5
subject <- rep(1:N, each = T)
NT <- N * T
Z1 <- rexp(NT, 1)
Z2 <- rbinom(NT, 1, 0.5)
Z <- cbind(Z1, Z2)
epsilon <- rnorm(NT, 0, 1)
X <- runif(NT, -2, 1)
psi <- c(-1, 0)
k <- length(psi)
PSI <- matrix(rep(psi, rep(NT, k)), ncol = k)
a <- rnorm(N, 0, 1)
A <- rep(a, each = T)
XP <- matrix(rep(X, k), nrow = NT)
XR <- cbind(1, X, pmax((XP - PSI), 0), Z)
bet <- c(1, -1, 3, -3, sqrt(3), -sqrt(3))
Y <- XR %*% bet + A + epsilon
tau <- 0.5
k <- 2

# Example 1: the working independence estimator; the error structure is "general"
est.type <- "WI"
wi.type <- "Compound"
tau <- 0.5
K.max <- 5
control <- fit.control(K.max = K.max)
Cn <- 1
mkqr.bea(y = Y, thre.x = X, cont.z = Z, id = subject, tau = tau, Cn = Cn,
         control = control, est.type = est.type, wi.type = wi.type)

# Example 2: the working correlated estimator; the correlation structure is "cs"
est.type <- "WC"
wc.type <- "cs"
tau <- 0.5
K.max <- 5
control <- fit.control(K.max = K.max)
Cn <- 1
mkqr.bea(y = Y, thre.x = X, cont.z = Z, id = subject, tau = tau, Cn = Cn,
         control = control, est.type = est.type, wc.type = wc.type)

## End(Not run)

Fit the multi-kink quantile regression conditional on a given or pre-specified number of change points.

Description

Fit the multi-kink quantile regression conditional on a given or pre-specified number of change points.

Usage

mkqr.fit(
  y,
  thre.x,
  cont.z,
  id,
  tau = 0.5,
  k,
  psi = NULL,
  bandwidth_type = "Hall-Sheather",
  control = fit.control(),
  est.type = "WI",
  wi.type = "general",
  wc.type = "cs"
)
mkqr.fit(
  y,
  thre.x,
  cont.z,
  id,
  tau = 0.5,
  k,
  psi = NULL,
  bandwidth_type = "Hall-Sheather",
  control = fit.control(),
  est.type = "WI",
  wi.type = "general",
  wc.type = "cs"
)

Arguments

`y`	A numeric vector of response.
`thre.x`	A numeric vector of scalar covariate with threshold effect.
`cont.z`	A numeric matrix of design with constant slopes.
`id`	A numeric vector of index used for longitudinal data; can be missing or NULL for iid data.
`tau`	The quantile level that belongs to (0,1).
`k`	The pre-specified number of change points.
`psi`	Numeric vector to indicate the starting values for the change points. When psi=NULL (default), k quantiles are assumed.
`bandwidth_type`	The bandwidth type. Specify one from "Hall-Sheather", "Bofinger", and "Chamberlain". Default is "Hall-Sheather".
`control`	A list returned by fit.control.
`est.type`	The estimation type for the longitudinal data. Specify one from "WI", "WC", corresponding to the working independence (WI) estimator and working correlation (WC) estimator. Default is "WI".
`wi.type`	If est.type = "WI", then set the error structure of the variance-covariance matrix estimation. Specify one from "Compound", "AR" and "general".
`wc.type`	If est.type = "WC", then set the correlation structure within subject. Specify one from "ar" and "cs". Default is "cs".

Value

A list containing the estimated regression coefficients with intercept (bet.est), the estimated standard error of the regression coefficients (bet.se), the estimated change points (psi.est), the estimated standard errors of threshold parameters (psi.se), and the fitted quantile objective value (rho).

Examples

## Not run: 
# Simple examples for iid data type
n <- 500
Z1 <- rexp(n, 1)
Z2 <- rbinom(n, 1, 0.5)
Z <- cbind(Z1, Z2)
epsilon <- rnorm(n, 0, 1)
X <- runif(n, -2, 1)
psi <- c(-1, 0)
k <- length(psi)
Y <- XR %*% bet + epsilon
result_iid <- mkqr.fit(Y, X, Z, tau = 0.5, k = k)

# Simple examples for longitudinal data
N <- 200
T <- 5
subject <- rep(1:N, each = T)
NT <- N * T
Z1 <- rexp(NT, 1)
Z2 <- rbinom(NT, 1, 0.5)
Z <- cbind(Z1, Z2)
epsilon <- rnorm(NT, 0, 1)
X <- runif(NT, -2, 1)
psi <- c(-1, 0)
k <- length(psi)
PSI <- matrix(rep(psi, rep(NT, k)), ncol = k)
a <- rnorm(N, 0, 1)
A <- rep(a, each = T)
XP <- matrix(rep(X, k), nrow = NT)
XR <- cbind(1, X, pmax((XP - PSI), 0), Z)
bet <- c(1, -1, 3, -3, sqrt(3), -sqrt(3))
Y <- XR %*% bet + A + epsilon
tau = 0.5
k = 2

# Example 1: the working independence estimator; the error structure is "general"
est.type = "WI";
wi.type = "Compound"
result_WI_Compound <- mkqr.fit(y = Y, thre.x = X, cont.z = Z, id = subject, tau = tau,
                               k = k, est.type = est.type, wi.type = wi.type)

# Example 2: the working correlated estimator; the correlation structure is "cs"
est.type = "WC";
wc.type = "cs"
result_WC_cs <- mkqr.fit(y = Y, thre.x = X, cont.z = Z, id = subject, tau = tau,
                         k = k, est.type = est.type, wc.type = wc.type)


## End(Not run)
## Not run: 
# Simple examples for iid data type
n <- 500
Z1 <- rexp(n, 1)
Z2 <- rbinom(n, 1, 0.5)
Z <- cbind(Z1, Z2)
epsilon <- rnorm(n, 0, 1)
X <- runif(n, -2, 1)
psi <- c(-1, 0)
k <- length(psi)
Y <- XR %*% bet + epsilon
result_iid <- mkqr.fit(Y, X, Z, tau = 0.5, k = k)

# Simple examples for longitudinal data
N <- 200
T <- 5
subject <- rep(1:N, each = T)
NT <- N * T
Z1 <- rexp(NT, 1)
Z2 <- rbinom(NT, 1, 0.5)
Z <- cbind(Z1, Z2)
epsilon <- rnorm(NT, 0, 1)
X <- runif(NT, -2, 1)
psi <- c(-1, 0)
k <- length(psi)
PSI <- matrix(rep(psi, rep(NT, k)), ncol = k)
a <- rnorm(N, 0, 1)
A <- rep(a, each = T)
XP <- matrix(rep(X, k), nrow = NT)
XR <- cbind(1, X, pmax((XP - PSI), 0), Z)
bet <- c(1, -1, 3, -3, sqrt(3), -sqrt(3))
Y <- XR %*% bet + A + epsilon
tau = 0.5
k = 2

# Example 1: the working independence estimator; the error structure is "general"
est.type = "WI";
wi.type = "Compound"
result_WI_Compound <- mkqr.fit(y = Y, thre.x = X, cont.z = Z, id = subject, tau = tau,
                               k = k, est.type = est.type, wi.type = wi.type)

# Example 2: the working correlated estimator; the correlation structure is "cs"
est.type = "WC";
wc.type = "cs"
result_WC_cs <- mkqr.fit(y = Y, thre.x = X, cont.z = Z, id = subject, tau = tau,
                         k = k, est.type = est.type, wc.type = wc.type)


## End(Not run)

Triceps Skinfold Thickness Dataset

Description

The data are derived from an anthropometric study of 892 females under 50 years in three Gambian villages in West Africa.

Usage

data(triceps)
data(triceps)

Format

A data frame with 892 observations on the following 3 variables:

age: Age of respondents.
lntriceps: Log of the triceps skinfold thickness.
triceps: Triceps skinfold thickness.

Source

Cole T.J., Green P.J. (1992). Smoothing reference centile curves: the LMS method and penalized likelihood. Statistics in medicine, 11(10): 1305-1319.

References

Cole T.J., Green P.J. (1992). Smoothing reference centile curves: the LMS method and penalized likelihood. Statistics in medicine, 11(10): 1305-1319.

Examples

data(triceps)
## maybe str(triceps) ...
data(triceps)
## maybe str(triceps) ...

Package 'MultiKink'

Help Index

Auxiliary parameters to control the model fitting

Description

Usage

Arguments

Value

Examples

Test the existence of kink effect in the multi-kink quantile regression

Description

Usage

Arguments

Value

Examples

Fit the multi-kink quantile regression in the absence of the number of change points.

Description

Usage

Arguments

Value

Examples

Fit the multi-kink quantile regression conditional on a given or pre-specified number of change points.

Description

Usage

Arguments

Value

Examples

Triceps Skinfold Thickness Dataset

Description

Usage

Format

Source

References

Examples