| Title: | Regularized Linear Models |
|---|---|
| Description: | Algorithms compute robust estimators for loss functions in the concave convex (CC) family by the iteratively reweighted convex optimization (IRCO), an extension of the iteratively reweighted least squares (IRLS). The IRCO reduces the weight of the observation that leads to a large loss; it also provides weights to help identify outliers. Applications include robust (penalized) generalized linear models and robust support vector machines. The package also contains penalized Poisson, negative binomial, zero-inflated Poisson, zero-inflated negative binomial regression models and robust models with non-convex loss functions. Wang et al. (2014) <doi:10.1002/sim.6314>, Wang et al. (2015) <doi:10.1002/bimj.201400143>, Wang et al. (2016) <doi:10.1177/0962280214530608>, Wang (2021) <doi:10.1007/s11749-021-00770-2>, Wang (2024) <doi:10.1111/anzs.12409>. |
| Authors: | Zhu Wang, with contributions from Achim Zeileis, Simon Jackman, Brian Ripley, and Patrick Breheny |
| Maintainer: | Zhu Wang <[email protected]> |
| License: | GPL-2 |
| Version: | 0.4-2.26 |
| Built: | 2024-09-26 06:22:42 UTC |
| Source: | CRAN |
conduct backward stepwise variable elimination for zero inflated count regression from zeroinfl function
be.zeroinfl(object, data, dist=c("poisson", "negbin", "geometric"), alpha=0.05, trace=FALSE)be.zeroinfl(object, data, dist=c("poisson", "negbin", "geometric"), alpha=0.05, trace=FALSE)
object |
an object from function zeroinfl |
data |
argument controlling formula processing
via |
dist |
one of the distributions in zeroinfl function |
alpha |
significance level of variable elimination |
trace |
logical value, if TRUE, print detailed calculation results |
conduct backward stepwise variable elimination for zero inflated count regression from zeroinfl function
an object of zeroinfl with all variables having p-values less than the significance level alpha
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
Generic function for extracting an estimator for the bread of sandwiches.
breadReg(x, which, ...)breadReg(x, which, ...)
x |
a fitted model object. |
which |
which penalty parameter(s)? |
... |
arguments passed to methods. |
A matrix containing an estimator for the penalized second derivative of log-likelihood function.
Typically, this should be an matrix corresponding
to parameters. The rows and columns should be named
as in coef or terms, respectively.
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) breadReg(fm_zinb, which=which.min(fm_zinb$bic))data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) breadReg(fm_zinb, which=which.min(fm_zinb$bic))
In a UK hospital, 135 expectant mothers were surveyed on the decision of breastfeeding their babies or not, along with two-level predictive factors
data(breastfeed)data(breastfeed)
Stephane Heritier, Eva Cantoni, Samuel Copt and Maria-Pia Victoria-Fese (2009). Robust Methods in Biostatistics, John Wiley & Sons
data(breastfeed) str(breastfeed)data(breastfeed) str(breastfeed)
Compute concave function values
compute_g(z, cfun, s, delta=0.0001)compute_g(z, cfun, s, delta=0.0001)
z |
vector nonnegative values from |
cfun |
integer from 1-8, concave function as in |
s |
a numeric value, see details in |
delta |
a positive small value, see details in |
Concave function values
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
u <- seq(0, 4, by=0.01) z <- u^2/2 ### this is dfun res <- compute_g(z, cfun=1, s=1) plot(z, res, ylab="Weight", type="l", lwd=2, main=expression(paste("hcave", " (", sigma, "=1)", )))u <- seq(0, 4, by=0.01) z <- u^2/2 ### this is dfun res <- compute_g(z, cfun=1, s=1) plot(z, res, ylab="Weight", type="l", lwd=2, main=expression(paste("hcave", " (", sigma, "=1)", )))
Weight value from concave function
compute_wt(z, weights, cfun, s, delta=0.0001)compute_wt(z, weights, cfun, s, delta=0.0001)
z |
vector nonnegative values from |
weights |
optional numeric vector of weights. |
cfun |
integer from 1-8, concave function as in |
s |
a numeric value, see details in |
delta |
a positive small value, see details in |
Weight value from concave function
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
u <- seq(0, 4, by=0.01) z <- u^2/2 ### this is dfun res <- compute_wt(z, cfun=1, s=1) plot(z, res, ylab="Weight", type="l", lwd=2, main=expression(paste("hcave", " (", sigma, "=1)", )))u <- seq(0, 4, by=0.01) z <- u^2/2 ### this is dfun res <- compute_wt(z, cfun=1, s=1) plot(z, res, ylab="Weight", type="l", lwd=2, main=expression(paste("hcave", " (", sigma, "=1)", )))
convert glm object to class glmreg, which then can be used for other purposes
conv2glmreg(object, family=c("poisson", "negbin"))conv2glmreg(object, family=c("poisson", "negbin"))
object |
an object of class glm |
family |
one of families in glm class |
an object of class glmreg
Zhu Wang <[email protected]>
convert zeroinfl object to class zipath, which then can be used to predict new data
conv2zipath(object, family=c("poisson", "negbin", "geometric"))conv2zipath(object, family=c("poisson", "negbin", "geometric"))
object |
an object of class zeroinfl |
family |
one of families in zeroinfl class |
an object of class zipath
Zhu Wang <[email protected]>
Does k-fold cross-validation for glmreg, produces a plot,
and returns cross-validated log-likelihood values for lambda
## S3 method for class 'formula' cv.glmreg(formula, data, weights, offset=NULL, contrasts=NULL, ...) ## S3 method for class 'matrix' cv.glmreg(x, y, weights, offset=NULL, ...) ## Default S3 method: cv.glmreg(x, ...) ## S3 method for class 'cv.glmreg' plot(x,se=TRUE,ylab=NULL, main=NULL, width=0.02, col="darkgrey", ...) ## S3 method for class 'cv.glmreg' predict(object, newx, ...) ## S3 method for class 'cv.glmreg' coef(object,which=object$lambda.which, ...)## S3 method for class 'formula' cv.glmreg(formula, data, weights, offset=NULL, contrasts=NULL, ...) ## S3 method for class 'matrix' cv.glmreg(x, y, weights, offset=NULL, ...) ## Default S3 method: cv.glmreg(x, ...) ## S3 method for class 'cv.glmreg' plot(x,se=TRUE,ylab=NULL, main=NULL, width=0.02, col="darkgrey", ...) ## S3 method for class 'cv.glmreg' predict(object, newx, ...) ## S3 method for class 'cv.glmreg' coef(object,which=object$lambda.which, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
x |
|
y |
response |
weights |
Observation weights; defaults to 1 per observation |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
contrasts |
the contrasts corresponding to |
object |
object of |
newx |
Matrix of values at which predictions are to be made. Not
used for |
which |
Indices of the penalty parameter |
se |
logical value, if |
ylab |
ylab on y-axis |
main |
title of plot |
width |
width of lines |
col |
color of standard error curve |
... |
Other arguments that can be passed to |
The function runs glmreg nfolds+1 times; the
first to compute the lambda sequence, and then to
compute the fit with each of the folds omitted. The error or the log-likelihood value is
accumulated, and the average value and standard deviation over the
folds is computed. Note that cv.glmreg can be used to search for
values for alpha: it is required to call cv.glmreg with a fixed vector foldid for different values of alpha.
an object of class "cv.glmreg" is returned, which is a
list with the ingredients of the cross-validation fit.
fit |
a fitted glmreg object for the full data. |
residmat |
matrix of log-likelihood values with row values for |
bic |
matrix of BIC values with row values for |
cv |
The mean cross-validated log-likelihood values - a vector of length
|
cv.error |
estimate of standard error of |
foldid |
an optional vector of values between 1 and |
lambda |
a vector of |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
glmreg and plot, predict, and coef methods for "cv.glmreg" object.
data("bioChemists", package = "pscl") fm_pois <- cv.glmreg(art ~ ., data = bioChemists, family = "poisson") title("Poisson Family",line=2.5) predict(fm_pois, newx=bioChemists[,-1])[1:4] coef(fm_pois)data("bioChemists", package = "pscl") fm_pois <- cv.glmreg(art ~ ., data = bioChemists, family = "poisson") title("Poisson Family",line=2.5) predict(fm_pois, newx=bioChemists[,-1])[1:4] coef(fm_pois)
Internal function to conduct k-fold cross-validation for glmreg, produces a plot,
and returns cross-validated log-likelihood values for lambda
cv.glmreg_fit(x, y, weights, offset, lambda=NULL, balance=TRUE, family=c("gaussian", "binomial", "poisson", "negbin"), type=c("loss", "error"), nfolds=10, foldid, plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)cv.glmreg_fit(x, y, weights, offset, lambda=NULL, balance=TRUE, family=c("gaussian", "binomial", "poisson", "negbin"), type=c("loss", "error"), nfolds=10, foldid, plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)
x |
|
y |
response |
weights |
Observation weights; defaults to 1 per observation |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
lambda |
Optional user-supplied lambda sequence; default is
|
balance |
for |
family |
response variable distribution |
type |
cross-validation criteria. For |
nfolds |
number of folds >=3, default is 10 |
foldid |
an optional vector of values between 1 and |
plot.it |
a logical value, to plot the estimated log-likelihood values if |
se |
a logical value, to plot with standard errors. |
parallel, n.cores
|
a logical value, parallel computing or not with the number of CPU cores to use. The cross-validation loop will attempt to send different CV folds off to different cores. |
trace |
a logical value, print progress of cross validation or not |
... |
Other arguments that can be passed to |
The function runs glmreg nfolds+1 times; the
first to compute the lambda sequence, and then to
compute the fit with each of the folds omitted. The error or the log-likelihood value is
accumulated, and the average value and standard deviation over the
folds is computed. Note that cv.glmreg can be used to search for
values for alpha: it is required to call cv.glmreg with a fixed vector foldid for different values of alpha.
an object of class "cv.glmreg" is returned, which is a
list with the ingredients of the cross-validation fit.
fit |
a fitted glmreg object for the full data. |
residmat |
matrix of log-likelihood values with row values for |
cv |
The mean cross-validated log-likelihood values - a vector of
|
cv.error |
estimate of standard error of |
foldid |
an optional vector of values between 1 and |
lambda |
a vector of |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
glmreg and plot, predict, and coef methods for "cv.glmreg" object.
Does k-fold cross-validation for glmregNB, produces a plot,
and returns cross-validated log-likelihood values for lambda
cv.glmregNB(formula, data, weights, offset=NULL, lambda=NULL, nfolds=10, foldid, plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)cv.glmregNB(formula, data, weights, offset=NULL, lambda=NULL, nfolds=10, foldid, plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)
formula |
symbolic description of the model |
data |
arguments controlling formula processing
via |
weights |
Observation weights; defaults to 1 per observation |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
lambda |
Optional user-supplied lambda sequence; default is
|
nfolds |
number of folds - default is 10. Although |
foldid |
an optional vector of values between 1 and |
plot.it |
a logical value, to plot the estimated log-likelihood values if |
se |
a logical value, to plot with standard errors. |
n.cores |
The number of CPU cores to use. The cross-validation loop will attempt to send different CV folds off to different cores. |
trace |
a logical value, print progress of cross-validation or not |
parallel |
a logical value, parallel computing or not |
... |
Other arguments that can be passed to |
The function runs glmregNB nfolds+1 times; the
first to get the lambda sequence, and then the remainder to
compute the fit with each of the folds omitted. The error is
accumulated, and the average error and standard deviation over the
folds is computed.
Note that cv.glmregNB does NOT search for
values for alpha. A specific value should be supplied, else
alpha=1 is assumed by default. If users would like to
cross-validate alpha as well, they should call cv.glmregNB
with a pre-computed vector foldid, and then use this same fold vector
in separate calls to cv.glmregNB with different values of
alpha.
an object of class "cv.glmregNB" is returned, which is a
list with the ingredients of the cross-validation fit.
fit |
a fitted glmregNB object for the full data. |
residmat |
matrix of log-likelihood values with row values for |
cv |
The mean cross-validated log-likelihood values - a vector of length
|
cv.error |
The standard error of cross-validated log-likelihood values - a vector of length
|
lambda |
a vector of |
foldid |
indicators of data used in each cross-validation, for reproductive purposes |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
glmregNB and plot, predict, and coef methods for "cv.glmregNB" object.
## Not run: data("bioChemists", package = "pscl") fm_nb <- cv.glmregNB(art ~ ., data = bioChemists) plot(fm_nb) ## End(Not run)## Not run: data("bioChemists", package = "pscl") fm_nb <- cv.glmregNB(art ~ ., data = bioChemists) plot(fm_nb) ## End(Not run)
Does k-fold cross-validation for irglmreg, produces a plot,
and returns cross-validated log-likelihood values for lambda
## S3 method for class 'formula' cv.irglmreg(formula, data, weights, offset=NULL, ...) ## S3 method for class 'matrix' cv.irglmreg(x, y, weights, offset=NULL, ...) ## Default S3 method: cv.irglmreg(x, ...) ## S3 method for class 'cv.irglmreg' plot(x,se=TRUE,ylab=NULL, main=NULL, width=0.02, col="darkgrey", ...) ## S3 method for class 'cv.irglmreg' coef(object,which=object$lambda.which, ...)## S3 method for class 'formula' cv.irglmreg(formula, data, weights, offset=NULL, ...) ## S3 method for class 'matrix' cv.irglmreg(x, y, weights, offset=NULL, ...) ## Default S3 method: cv.irglmreg(x, ...) ## S3 method for class 'cv.irglmreg' plot(x,se=TRUE,ylab=NULL, main=NULL, width=0.02, col="darkgrey", ...) ## S3 method for class 'cv.irglmreg' coef(object,which=object$lambda.which, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
x |
|
y |
response |
weights |
Observation weights; defaults to 1 per observation |
offset |
Not implemented yet |
object |
object of |
which |
Indices of the penalty parameter |
se |
logical value, if |
ylab |
ylab on y-axis |
main |
title of plot |
width |
width of lines |
col |
color of standard error curve |
... |
Other arguments that can be passed to |
The function runs irglmreg nfolds+1 times; the
first to compute the lambda sequence, and then to
compute the fit with each of the folds omitted. The error or the loss value is
accumulated, and the average value and standard deviation over the
folds is computed. Note that cv.irglmreg can be used to search for
values for alpha: it is required to call cv.irglmreg with a fixed vector foldid for different values of alpha.
an object of class "cv.irglmreg" is returned, which is a
list with the ingredients of the cross-validation fit.
fit |
a fitted irglmreg object for the full data. |
residmat |
matrix of log-likelihood values with row values for |
bic |
matrix of BIC values with row values for |
cv |
The mean cross-validated log-likelihood values - a vector of length
|
cv.error |
estimate of standard error of |
foldid |
an optional vector of values between 1 and |
lambda |
a vector of |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
irglmreg and plot, predict, and coef methods for "cv.irglmreg" object.
Internal function to conduct k-fold cross-validation for irglmreg, produces a plot,
and returns cross-validated loss values for lambda
cv.irglmreg_fit(x, y, weights, offset, lambda=NULL, balance=TRUE, cfun=4, dfun=1, s=1.5, nfolds=10, foldid, type = c("loss", "error"), plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)cv.irglmreg_fit(x, y, weights, offset, lambda=NULL, balance=TRUE, cfun=4, dfun=1, s=1.5, nfolds=10, foldid, type = c("loss", "error"), plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)
x |
|
y |
response |
weights |
Observation weights; defaults to 1 per observation |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
lambda |
Optional user-supplied lambda sequence; default is
|
balance |
for |
cfun |
a number from 1 to 7, type of convex cap (concave) function |
dfun |
a number from 1, 4-7, type of convex downward function |
s |
nonconvex loss tuning parameter for robust regression and classification. |
nfolds |
number of folds >=3, default is 10 |
foldid |
an optional vector of values between 1 and |
type |
cross-validation criteria. For |
plot.it |
a logical value, to plot the estimated log-likelihood values if |
se |
a logical value, to plot with standard errors. |
n.cores |
The number of CPU cores to use. The cross-validation loop will attempt to send different CV folds off to different cores. |
trace |
a logical value, print progress of cross validation or not |
parallel |
a logical value, parallel computing or not |
... |
Other arguments that can be passed to |
The function runs irglmreg nfolds+1 times; the
first to compute the lambda sequence, and then to
compute the fit with each of the folds omitted. The error or the log-likelihood value is
accumulated, and the average value and standard deviation over the
folds is computed. Note that cv.irglmreg can be used to search for
values for alpha: it is required to call cv.irglmreg with a fixed vector foldid for different values of alpha.
an object of class "cv.irglmreg" is returned, which is a
list with the ingredients of the cross-validation fit.
fit |
a fitted irglmreg object for the full data. |
residmat |
matrix of loss values or errors with row values for |
cv |
The mean cross-validated loss values or errors - a vector of length
|
cv.error |
estimate of standard error of |
foldid |
an optional vector of values between 1 and |
lambda |
a vector of |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
irglmreg and plot, predict, and coef methods for "cv.irglmreg" object.
Does k-fold cross-validation for irsvm
## S3 method for class 'formula' cv.irsvm(formula, data, weights, contrasts=NULL, ...) ## S3 method for class 'matrix' cv.irsvm(x, y, weights, ...) ## Default S3 method: cv.irsvm(x, ...)## S3 method for class 'formula' cv.irsvm(formula, data, weights, contrasts=NULL, ...) ## S3 method for class 'matrix' cv.irsvm(x, y, weights, ...) ## Default S3 method: cv.irsvm(x, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
x |
|
y |
response |
weights |
Observation weights; defaults to 1 per observation |
contrasts |
the contrasts corresponding to |
... |
Other arguments that can be passed to |
Does a K-fold cross-validation to determine optimal tuning parameters in SVM: cost and gamma if kernel is nonlinear. It can also choose s used in cfun.
An object contains a list of ingredients of cross-validation including optimal tuning parameters.
residmat |
matrix with row values for |
cost |
a value of |
gamma |
a value of |
s |
value of |
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
## Not run: x <- matrix(rnorm(40*2), ncol=2) y <- c(rep(-1, 20), rep(1, 20)) x[y==1,] <- x[y==1, ] + 1 irsvm.opt <- cv.irsvm(x, y, type="C-classification", s=1, kernel="linear", cfun="acave") irsvm.opt$cost irsvm.opt$gamma irsvm.opt$s ## End(Not run)## Not run: x <- matrix(rnorm(40*2), ncol=2) y <- c(rep(-1, 20), rep(1, 20)) x[y==1,] <- x[y==1, ] + 1 irsvm.opt <- cv.irsvm(x, y, type="C-classification", s=1, kernel="linear", cfun="acave") irsvm.opt$cost irsvm.opt$gamma irsvm.opt$s ## End(Not run)
Internal function to conduct k-fold cross-validation for irsvm
cv.irsvm_fit(x, y, weights, cfun="ccave", s=c(1, 5), type=NULL, kernel="radial", gamma=2^(-4:10), cost=2^(-4:4), epsilon=0.1, balance=TRUE, nfolds=10, foldid, trim_ratio=0.9, n.cores=2, ...)cv.irsvm_fit(x, y, weights, cfun="ccave", s=c(1, 5), type=NULL, kernel="radial", gamma=2^(-4:10), cost=2^(-4:4), epsilon=0.1, balance=TRUE, nfolds=10, foldid, trim_ratio=0.9, n.cores=2, ...)
x |
a data matrix, a vector, or a sparse 'design matrix' (object of class
|
y |
a response vector with one label for each row/component of
|
weights |
the weight of each subject. It should be in the same length of |
cfun |
character, type of convex cap (concave) function.
|
s |
tuning parameter of |
type |
|
kernel, gamma
|
the kernel used in training and predicting. You
might consider changing some of the following parameters, depending
on the kernel type.
|
cost |
cost of constraints violation (default: 1)—it is the
‘C’-constant of the regularization term in the Lagrange formulation. This is proportional to the inverse of |
epsilon |
epsilon in the insensitive-loss function (default: 0.1) |
balance |
for |
nfolds |
number of folds >=3, default is 10 |
foldid |
an optional vector of values between 1 and |
trim_ratio |
a number between 0 and 1 for trimmed least squares, useful if |
n.cores |
The number of CPU cores to use. The cross-validation loop will attempt to send different CV folds off to different cores. |
... |
Other arguments that can be passed to |
This function is the driving force behind cv.irsvm. Does a K-fold cross-validation to determine optimal tuning parameters in SVM: cost and gamma if kernel is nonlinear. It can also choose s used in cfun.
an object of class "cv.irsvm" is returned, which is a
list with the ingredients of the cross-validation fit.
residmat |
matrix with row values for |
cost |
a value of |
gamma |
a value of |
s |
value of |
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
Does k-fold cross-validation for nclreg, produces a plot,
and returns cross-validated loss values for lambda
## S3 method for class 'formula' cv.nclreg(formula, data, weights, offset=NULL, ...) ## S3 method for class 'matrix' cv.nclreg(x, y, weights, offset=NULL, ...) ## Default S3 method: cv.nclreg(x, ...) ## S3 method for class 'cv.nclreg' plot(x,se=TRUE,ylab=NULL, main=NULL, width=0.02, col="darkgrey", ...) ## S3 method for class 'cv.nclreg' coef(object,which=object$lambda.which, ...)## S3 method for class 'formula' cv.nclreg(formula, data, weights, offset=NULL, ...) ## S3 method for class 'matrix' cv.nclreg(x, y, weights, offset=NULL, ...) ## Default S3 method: cv.nclreg(x, ...) ## S3 method for class 'cv.nclreg' plot(x,se=TRUE,ylab=NULL, main=NULL, width=0.02, col="darkgrey", ...) ## S3 method for class 'cv.nclreg' coef(object,which=object$lambda.which, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
x |
|
y |
response |
weights |
Observation weights; defaults to 1 per observation |
offset |
Not implemented yet |
object |
object of |
which |
Indices of the penalty parameter |
se |
logical value, if |
ylab |
ylab on y-axis |
main |
title of plot |
width |
width of lines |
col |
color of standard error curve |
... |
Other arguments that can be passed to |
The function runs nclreg nfolds+1 times; the
first to compute the lambda sequence, and then to
compute the fit with each of the folds omitted. The error or the loss value is
accumulated, and the average value and standard deviation over the
folds is computed. Note that cv.nclreg can be used to search for
values for alpha: it is required to call cv.nclreg with a fixed vector foldid for different values of alpha.
an object of class "cv.nclreg" is returned, which is a
list with the ingredients of the cross-validation fit.
fit |
a fitted nclreg object for the full data. |
residmat |
matrix of loss values with row values for |
bic |
matrix of BIC values with row values for |
cv |
The mean cross-validated loss values - a vector of length
|
cv.error |
estimate of standard error of |
foldid |
an optional vector of values between 1 and |
lambda |
a vector of |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang (2021), MM for Penalized Estimation, TEST, doi:10.1007/s11749-021-00770-2
nclreg and plot, predict, and coef methods for "cv.nclreg" object.
Internal function to conduct k-fold cross-validation for nclreg, produces a plot,
and returns cross-validated loss values for lambda
cv.nclreg_fit(x, y, weights, offset, lambda=NULL, balance=TRUE, rfamily=c("clossR", "closs", "gloss", "qloss"), s=1.5, nfolds=10, foldid, type = c("loss", "error"), plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)cv.nclreg_fit(x, y, weights, offset, lambda=NULL, balance=TRUE, rfamily=c("clossR", "closs", "gloss", "qloss"), s=1.5, nfolds=10, foldid, type = c("loss", "error"), plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)
x |
|
y |
response |
weights |
Observation weights; defaults to 1 per observation |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
lambda |
Optional user-supplied lambda sequence; default is
|
balance |
for |
rfamily |
response variable distribution and nonconvex loss function |
s |
nonconvex loss tuning parameter for robust regression and classification. |
nfolds |
number of folds >=3, default is 10 |
foldid |
an optional vector of values between 1 and |
type |
cross-validation criteria. For |
plot.it |
a logical value, to plot the estimated loss values if |
se |
a logical value, to plot with standard errors. |
n.cores |
The number of CPU cores to use. The cross-validation loop will attempt to send different CV folds off to different cores. |
trace |
a logical value, print progress of cross validation or not |
parallel |
a logical value, parallel computing or not |
... |
Other arguments that can be passed to |
The function runs nclreg nfolds+1 times; the
first to compute the lambda sequence, and then to
compute the fit with each of the folds omitted. The error or the loss value is
accumulated, and the average value and standard deviation over the
folds is computed. Note that cv.nclreg can be used to search for
values for alpha: it is required to call cv.nclreg with a fixed vector foldid for different values of alpha.
an object of class "cv.nclreg" is returned, which is a
list with the ingredients of the cross-validation fit.
fit |
a fitted nclreg object for the full data. |
residmat |
matrix of loss values or errors with row values for |
cv |
The mean cross-validated loss values or errors - a vector of length
|
cv.error |
estimate of standard error of |
foldid |
an optional vector of values between 1 and |
lambda |
a vector of |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang (2021), MM for Penalized Estimation, TEST, doi:10.1007/s11749-021-00770-2
nclreg and plot, predict, and coef methods for "cv.nclreg" object.
Does k-fold cross-validation for zipath, produces a plot,
and returns cross-validated log-likelihood values for lambda
## S3 method for class 'formula' cv.zipath(formula, data, weights, offset=NULL, contrasts=NULL, ...) ## S3 method for class 'matrix' cv.zipath(X, Z, Y, weights, offsetx=NULL, offsetz=NULL, ...) ## Default S3 method: cv.zipath(X, ...) ## S3 method for class 'cv.zipath' predict(object, newdata, ...) ## S3 method for class 'cv.zipath' coef(object, which=object$lambda.which, model = c("full", "count", "zero"), ...)## S3 method for class 'formula' cv.zipath(formula, data, weights, offset=NULL, contrasts=NULL, ...) ## S3 method for class 'matrix' cv.zipath(X, Z, Y, weights, offsetx=NULL, offsetz=NULL, ...) ## Default S3 method: cv.zipath(X, ...) ## S3 method for class 'cv.zipath' predict(object, newdata, ...) ## S3 method for class 'cv.zipath' coef(object, which=object$lambda.which, model = c("full", "count", "zero"), ...)
formula |
symbolic description of the model with an optional numeric vector |
data |
arguments controlling formula processing
via |
weights |
Observation weights; defaults to 1 per observation |
offset |
optional numeric vector with an a priori known component to be included in the linear predictor of the count model or zero model. See below for an example. |
X |
predictor matrix of the count model |
Z |
predictor matrix of the zero model |
Y |
response variable |
offsetx, offsetz
|
optional numeric vector with an a priori known component to be included in the linear predictor of the count model (offsetx)or zero model (offsetz). |
contrasts |
a list with elements |
object |
object of class |
newdata |
optionally, a data frame in which to look for variables with which to predict. If omitted, the original observations are used. |
which |
Indices of the pair of penalty parameters |
model |
character specifying for which component of the model the estimated coefficients should be extracted. |
... |
Other arguments that can be passed to |
The function runs zipath nfolds+1 times; the
first to compute the (lambda.count, lambda.zero) sequence, and then to
compute the fit with each of the folds omitted.
The model is fitted to the training data and then given the fitted model the log-likelihood is evaluated at the observations left out, i.e., the test data.
The average value of log-likelihood and standard deviation over the
folds is computed. Note that cv.zipath can be used to search for
values for count.alpha or zero.alpha: it is required to call cv.zipath with a fixed vector foldid for different values of count.alpha or zero.alpha.
The methods for coef and predict were deprecated since version 0.3-25. In fact, the fit object was removed in the output of cv.zipath so that predict an object of cv.zipath is not feasible, and should be via zipath. See examples below. The reason for such a change is that cv.zipath can take both formula and matrix, hence predict on cv. zipath object can easily lead to problems in codes.
When family="negbin", it can be slow because there is a repeated search for the theta values by default. One may change the default values from init.theta=NULL, theta.fixed=FALSE to init.theta=MLE, theta.fixed=TRUE, where MLE is a number from glm.nb in the R package MASS or something desired.
an object of class "cv.zipath" is returned, which is a
list with the components of the cross-validation fit.
fit |
a fitted zipath object for the full data. |
residmat |
matrix for cross-validated log-likelihood at each |
bic |
matrix of BIC values with row values for |
cv |
The mean cross-validated log-likelihood - a vector of length
|
cv.error |
estimate of standard error of |
foldid |
an optional vector of values between 1 and |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
zipath and plot, predict, methods for "cv.zipath" object.
## Not run: data("bioChemists", package = "pscl") fm_zip <- zipath(art ~ . | ., data = bioChemists, family = "poisson", nlambda=10) fm_cvzip <- cv.zipath(art ~ . | ., data = bioChemists, family = "poisson", nlambda=10) ### prediction from the best model pred <- predict(fm_zip, newdata=bioChemists, which=fm_cvzip$lambda.which) coef(fm_zip, which=fm_cvzip$lambda.which) fm_znb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) fm_cvznb <- cv.zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) pred <- predict(fm_znb, which=fm_cvznb$lambda.which) coef(fm_znb, which=fm_cvznb$lambda.which) fm_zinb2 <- zipath(art ~ . +offset(log(phd))| ., data = bioChemists, family = "poisson", nlambda=10) fm_cvzinb2 <- cv.zipath(art ~ . +offset(log(phd))| ., data = bioChemists, family = "poisson", nlambda=10) pred <- predict(fm_zinb2, which=fm_cvzinb2$lambda.which) coef(fm_zinb2, which=fm_cvzinb2$lambda.which) ## End(Not run)## Not run: data("bioChemists", package = "pscl") fm_zip <- zipath(art ~ . | ., data = bioChemists, family = "poisson", nlambda=10) fm_cvzip <- cv.zipath(art ~ . | ., data = bioChemists, family = "poisson", nlambda=10) ### prediction from the best model pred <- predict(fm_zip, newdata=bioChemists, which=fm_cvzip$lambda.which) coef(fm_zip, which=fm_cvzip$lambda.which) fm_znb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) fm_cvznb <- cv.zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) pred <- predict(fm_znb, which=fm_cvznb$lambda.which) coef(fm_znb, which=fm_cvznb$lambda.which) fm_zinb2 <- zipath(art ~ . +offset(log(phd))| ., data = bioChemists, family = "poisson", nlambda=10) fm_cvzinb2 <- cv.zipath(art ~ . +offset(log(phd))| ., data = bioChemists, family = "poisson", nlambda=10) pred <- predict(fm_zinb2, which=fm_cvzinb2$lambda.which) coef(fm_zinb2, which=fm_cvzinb2$lambda.which) ## End(Not run)
Internal function k-fold cross-validation for zipath, produces a plot, and returns cross-validated log-likelihood values for lambda
cv.zipath_fit(X, Z, Y, weights, offsetx, offsetz, nlambda=100, lambda.count=NULL, lambda.zero=NULL, nfolds=10, foldid, plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)cv.zipath_fit(X, Z, Y, weights, offsetx, offsetz, nlambda=100, lambda.count=NULL, lambda.zero=NULL, nfolds=10, foldid, plot.it=TRUE, se=TRUE, n.cores=2, trace=FALSE, parallel=FALSE, ...)
X |
predictor matrix of the count model |
Z |
predictor matrix of the zero model |
Y |
response variable |
weights |
optional numeric vector of weights. |
offsetx |
optional numeric vector with an a priori known component to be included in the linear predictor of the count model. |
offsetz |
optional numeric vector with an a priori known component to be included in the linear predictor of the zero model. |
nlambda |
number of |
lambda.count |
Optional user-supplied lambda.count sequence; default is
|
lambda.zero |
Optional user-supplied lambda.zero sequence; default is
|
nfolds |
number of folds >=3, default is 10 |
foldid |
an optional vector of values between 1 and |
plot.it |
a logical value, to plot the estimated log-likelihood values if |
se |
a logical value, to plot with standard errors. |
n.cores |
The number of CPU cores to use. The cross-validation loop will attempt to send different CV folds off to different cores. |
trace |
a logical value, print progress of cross-validation or not |
parallel |
a logical value, parallel computing or not |
... |
Other arguments that can be passed to |
The function runs zipath nfolds+1 times; the
first to compute the (lambda.count, lambda.zero) sequence, and then to
compute the fit with each of the folds omitted. The log-likelihood value is
accumulated, and the average value and standard deviation over the
folds is computed. Note that cv.zipath can be used to search for
values for count.alpha or zero.alpha: it is required to call cv.zipath with a fixed vector foldid for different values of count.alpha or zero.alpha.
The method for coef by default
return a single vector of coefficients, i.e., all coefficients are concatenated. By setting the model
argument, the estimates for the corresponding model components can be extracted.
an object of class "cv.zipath" is returned, which is a
list with the components of the cross-validation fit.
fit |
a fitted zipath object for the full data. |
residmat |
matrix for cross-validated log-likelihood at each |
bic |
matrix of BIC values with row values for |
cv |
The mean cross-validated log-likelihood - a vector of length
|
cv.error |
estimate of standard error of |
foldid |
an optional vector of values between 1 and |
lambda.which |
index of |
lambda.optim |
value of |
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
zipath and plot, predict, and coef methods for "cv.zipath" object.
A cohort of 3066 Americans over the age of 50 were studied on health care utilization, doctor office visits.
data(docvisits)data(docvisits)
Stephane Heritier, Eva Cantoni, Samuel Copt and Maria-Pia Victoria-Fese (2009). Robust Methods in Biostatistics, John Wiley & Sons
data(docvisits) str(docvisits)data(docvisits) str(docvisits)
Generic function for extracting the empirical first derivative of log-likelihood function of a fitted regularized model.
estfunReg(x, ...)estfunReg(x, ...)
x |
a fitted model object. |
... |
arguments passed to methods. |
A matrix containing the empirical first derivative of log-likelihood functions.
Typically, this should be an matrix corresponding
to observations and parameters. The columns should be named
as in coef or terms, respectively.
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) res <- estfunReg(fm_zinb, which=which.min(fm_zinb$bic))data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) res <- estfunReg(fm_zinb, which=which.min(fm_zinb$bic))
Compute response value to raw prediction such as linear predictor in GLM
gfunc(mu, family, epsbino)gfunc(mu, family, epsbino)
mu |
vector of numbers as response value in GLM, for instance, probability estimation if |
family |
integer from 1-4, corresponding to "gaussian", "binomial", "poisson", "negbin", respectively |
epsbino |
a small positive value for |
linear predictor f=x'b for predictor x and coefficient b if the model is linear
Fit a generalized linear model via penalized maximum likelihood. The regularization path is computed for the lasso (or elastic net penalty), scad (or snet) and mcp (or mnet penalty), at a grid of values for the regularization parameter lambda. Fits linear, logistic, Poisson and negative binomial (fixed scale parameter) regression models.
## S3 method for class 'formula' glmreg(formula, data, weights, offset=NULL, contrasts=NULL, x.keep=FALSE, y.keep=TRUE, ...) ## S3 method for class 'matrix' glmreg(x, y, weights, offset=NULL, ...) ## Default S3 method: glmreg(x, ...)## S3 method for class 'formula' glmreg(formula, data, weights, offset=NULL, contrasts=NULL, x.keep=FALSE, y.keep=TRUE, ...) ## S3 method for class 'matrix' glmreg(x, y, weights, offset=NULL, ...) ## Default S3 method: glmreg(x, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights. If |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
x |
input matrix, of dimension nobs x nvars; each row is an observation vector |
y |
response variable. Quantitative for |
x.keep, y.keep
|
logical values: keep response variables or keep response variable? |
contrasts |
the contrasts corresponding to |
... |
Other arguments passing to |
The sequence of models implied by lambda is fit by coordinate
descent. For family="gaussian" this is the lasso, mcp or scad sequence if
alpha=1, else it is the enet, mnet or snet sequence.
For the other families, this is a lasso (mcp, scad) or elastic net (mnet, snet) regularization path
for fitting the generalized linear regression
paths, by maximizing the appropriate penalized log-likelihood.
Note that the objective function for "gaussian" is
if standardize=FALSE and
if standardize=TRUE. For the other models it is
if standardize=FALSE and
if standardize=TRUE.
An object with S3 class "glmreg" for the various types of models.
call |
the call that produced this object |
b0 |
Intercept sequence of length |
beta |
A |
lambda |
The actual sequence of |
offset |
the offset vector used. |
resdev |
The computed deviance (for |
nulldev |
Null deviance (per observation). This is defined to be 2*(loglike_sat -loglike(Null)); The NULL model refers to the intercept model. |
nobs |
number of observations |
pll |
penalized log-likelihood values for standardized coefficients in the IRLS iterations. For |
pllres |
penalized log-likelihood value for the estimated model on the original scale of coefficients |
fitted.values |
the fitted mean values, obtained by transforming the linear predictors by the inverse of the link function. |
Zhu Wang <[email protected]>
Breheny, P. and Huang, J. (2011) Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection. Ann. Appl. Statist., 5: 232-253.
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
print, predict, coef and plot methods, and the cv.glmreg function.
#binomial x=matrix(rnorm(100*20),100,20) g2=sample(0:1,100,replace=TRUE) fit2=glmreg(x,g2,family="binomial") #poisson and negative binomial data("bioChemists", package = "pscl") fm_pois <- glmreg(art ~ ., data = bioChemists, family = "poisson") coef(fm_pois) fm_nb1 <- glmreg(art ~ ., data = bioChemists, family = "negbin", theta=1) coef(fm_nb1) #offset x <- matrix(rnorm(100*20),100,20) y <- rpois(100, lambda=1) exposure <- rep(0.5, length(y)) fit2 <- glmreg(x,y, lambda=NULL, nlambda=10, lambda.min.ratio=1e-4, offset=log(exposure), family="poisson") predict(fit2, newx=x, newoffset=log(exposure)) ## Not run: fm_nb2 <- glmregNB(art ~ ., data = bioChemists) coef(fm_nb2) ## End(Not run)#binomial x=matrix(rnorm(100*20),100,20) g2=sample(0:1,100,replace=TRUE) fit2=glmreg(x,g2,family="binomial") #poisson and negative binomial data("bioChemists", package = "pscl") fm_pois <- glmreg(art ~ ., data = bioChemists, family = "poisson") coef(fm_pois) fm_nb1 <- glmreg(art ~ ., data = bioChemists, family = "negbin", theta=1) coef(fm_nb1) #offset x <- matrix(rnorm(100*20),100,20) y <- rpois(100, lambda=1) exposure <- rep(0.5, length(y)) fit2 <- glmreg(x,y, lambda=NULL, nlambda=10, lambda.min.ratio=1e-4, offset=log(exposure), family="poisson") predict(fit2, newx=x, newoffset=log(exposure)) ## Not run: fm_nb2 <- glmregNB(art ~ ., data = bioChemists) coef(fm_nb2) ## End(Not run)
Fit a generalized linear model via penalized maximum likelihood. The regularization path is computed for the lasso (or elastic net penalty), snet and mnet penalty, at a grid of values for the regularization parameter lambda. Fits linear, logistic, Poisson and negative binomial (fixed scale parameter) regression models.
glmreg_fit(x, y, weights, start=NULL, etastart=NULL, mustart=NULL, offset = NULL, nlambda=100, lambda=NULL, lambda.min.ratio=ifelse(nobs<nvars,.05, .001), alpha=1, gamma=3, rescale=TRUE, standardize=TRUE, intercept=TRUE, penalty.factor = rep(1, nvars), thresh=1e-6, eps.bino=1e-5, maxit=1000, eps=.Machine$double.eps, theta, family=c("gaussian", "binomial", "poisson", "negbin"), penalty=c("enet","mnet","snet"), convex=FALSE, x.keep=FALSE, y.keep=TRUE, trace=FALSE)glmreg_fit(x, y, weights, start=NULL, etastart=NULL, mustart=NULL, offset = NULL, nlambda=100, lambda=NULL, lambda.min.ratio=ifelse(nobs<nvars,.05, .001), alpha=1, gamma=3, rescale=TRUE, standardize=TRUE, intercept=TRUE, penalty.factor = rep(1, nvars), thresh=1e-6, eps.bino=1e-5, maxit=1000, eps=.Machine$double.eps, theta, family=c("gaussian", "binomial", "poisson", "negbin"), penalty=c("enet","mnet","snet"), convex=FALSE, x.keep=FALSE, y.keep=TRUE, trace=FALSE)
x |
input matrix, of dimension nobs x nvars; each row is an observation vector. |
y |
response variable. Quantitative for |
weights |
observation weights. Can be total counts if responses are proportion matrices. Default is 1 for each observation |
start |
starting values for the parameters in the linear predictor. |
etastart |
starting values for the linear predictor. |
mustart |
starting values for the vector of means. |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
nlambda |
The number of |
lambda |
by default, the algorithm provides a sequence of regularization values, or a user supplied |
lambda.min.ratio |
Smallest value for |
alpha |
The |
gamma |
The tuning parameter of the |
rescale |
logical value, if TRUE, adaptive rescaling of the penalty parameter for |
standardize |
logical value for x variable standardization, prior to
fitting the model sequence. The coefficients are always returned on
the original scale. Default is |
intercept |
logical value: if TRUE (default), intercept(s) are fitted; otherwise, intercept(s) are set to zero |
penalty.factor |
This is a number that multiplies |
thresh |
Convergence threshold for coordinate descent. Defaults value is |
eps.bino |
a lower bound of probabilities to be truncated, for computing weights and related values when |
maxit |
Maximum number of coordinate descent iterations for each |
eps |
If a coefficient is less than |
convex |
Calculate index for which objective function ceases to
be locally convex? Default is FALSE and only useful if |
theta |
an overdispersion scaling parameter for |
family |
Response type (see above) |
penalty |
Type of regularization |
x.keep, y.keep
|
For glmreg: logical values indicating whether the response vector and model matrix used in the fitting process should be returned as components of the returned value. For glmreg_fit: x is a design matrix of dimension n * p, and x is a vector of observations of length n. |
trace |
If |
The sequence of models implied by lambda is fit by coordinate
descent. For family="gaussian" this is the lasso, mcp or scad sequence if
alpha=1, else it is the enet, mnet or snet sequence.
For the other families, this is a lasso (mcp, scad) or elastic net (mnet, snet) regularization path
for fitting the generalized linear regression
paths, by maximizing the appropriate penalized log-likelihood.
Note that the objective function for "gaussian" is
if standardize=FALSE and
if standardize=TRUE. For the other models it is
if standardize=FALSE and
if standardize=TRUE.
An object with S3 class "glmreg" for the various types of models.
call |
the call that produced the model fit |
b0 |
Intercept sequence of length |
beta |
A |
lambda |
The actual sequence of |
satu |
satu=1 if a saturated model (deviance/null deviance < 0.05) is fit. Otherwise satu=0. The number of |
dev |
The computed deviance (for |
nulldev |
Null deviance (per observation). This is defined to be 2*(loglike_sat -loglike(Null)); The NULL model refers to the intercept model. |
nobs |
number of observations |
Zhu Wang <[email protected]>
Breheny, P. and Huang, J. (2011) Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection. Ann. Appl. Statist., 5: 232-253.
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Fit a negative binomial linear model via penalized maximum likelihood. The regularization path is computed for the lasso (or elastic net penalty), snet and mnet penalty, at a grid of values for the regularization parameter lambda.
glmregNB(formula, data, weights, offset=NULL, nlambda = 100, lambda=NULL, lambda.min.ratio = ifelse(nobs<nvars,0.05,0.001), alpha=1, gamma=3, rescale=TRUE, standardize = TRUE, penalty.factor = rep(1, nvars), thresh = 0.001, maxit.theta = 10, maxit=1000, eps=.Machine$double.eps, trace=FALSE, start = NULL, etastart = NULL, mustart = NULL, theta.fixed=FALSE, theta0=NULL, init.theta=NULL, link=log, penalty=c("enet","mnet","snet"), method="glmreg_fit", model=TRUE, x.keep=FALSE, y.keep=TRUE, contrasts=NULL, convex=FALSE, parallel=TRUE, n.cores=2, ...)glmregNB(formula, data, weights, offset=NULL, nlambda = 100, lambda=NULL, lambda.min.ratio = ifelse(nobs<nvars,0.05,0.001), alpha=1, gamma=3, rescale=TRUE, standardize = TRUE, penalty.factor = rep(1, nvars), thresh = 0.001, maxit.theta = 10, maxit=1000, eps=.Machine$double.eps, trace=FALSE, start = NULL, etastart = NULL, mustart = NULL, theta.fixed=FALSE, theta0=NULL, init.theta=NULL, link=log, penalty=c("enet","mnet","snet"), method="glmreg_fit", model=TRUE, x.keep=FALSE, y.keep=TRUE, contrasts=NULL, convex=FALSE, parallel=TRUE, n.cores=2, ...)
formula |
formula used to describe a model. |
data |
argument controlling formula processing
via |
weights |
an optional vector of ‘prior weights’ to be used in the fitting process. Should be |
offset |
optional numeric vector with an a priori known component to be included in the linear predictor of the model. |
nlambda |
The number of |
lambda |
A user supplied |
lambda.min.ratio |
Smallest value for |
alpha |
The L2 penalty mixing parameter, with
|
gamma |
The tuning parameter of the |
rescale |
logical value, if TRUE, adaptive rescaling of the penalty parameter for |
standardize |
Logical flag for x variable standardization, prior to
fitting the model sequence. The coefficients are always returned on
the original scale. Default is |
penalty.factor |
This is a number that multiplies |
thresh |
Convergence threshold for coordinate descent. Defaults value is |
maxit.theta |
Maximum number of iterations for estimating |
maxit |
Maximum number of coordinate descent iterations for each |
eps |
If a number is less than |
trace |
If |
start, etastart, mustart, ...
|
arguments for the |
init.theta |
initial scaling parameter |
theta.fixed |
Estimate scale parameter theta? Default is FALSE. Note, the algorithm may become slow. In this case, one may use |
.
theta0 |
initial scale parameter vector theta, with length |
convex |
Calculate index for which objective function ceases to
be locally convex? Default is FALSE and only useful if |
link |
link function, default is |
penalty |
Type of regularization |
method |
estimation method |
model, x.keep, y.keep
|
logicals. If |
contrasts |
the contrasts corresponding to |
parallel, n.cores
|
a logical value, parallel computing or not for sequence of |
The sequence of models implied by lambda is fit by coordinate
descent. This is a lasso (mcp, scad) or elastic net (mnet, snet) regularization path
for fitting the negative binomial linear regression
paths, by maximizing the penalized log-likelihood.
Note that the objective function is
if standardize=FALSE and
if standardize=TRUE.
An object with S3 class "glmreg", "glmregNB" for the various types of models.
call |
the call that produced the model fit |
b0 |
Intercept sequence of length |
beta |
A |
lambda |
The actual sequence of |
resdev |
The computed deviance. The deviance calculations incorporate weights if present in the model. The deviance is defined to be 2*(loglike_sat - loglike), where loglike_sat is the log-likelihood for the saturated model (a model with a free parameter per observation). |
nulldev |
Null deviance (per observation). This is defined to be 2*(loglike_sat -loglike(Null)); The NULL model refers to the intercept model. |
nobs |
number of observations |
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
print, predict, coef and plot methods, and the cv.glmregNB function.
## Not run: data("bioChemists", package = "pscl") system.time(fm_nb1 <- glmregNB(art ~ ., data = bioChemists, parallel=FALSE)) system.time(fm_nb2 <- glmregNB(art ~ ., data = bioChemists, parallel=TRUE, n.cores=2)) coef(fm_nb1) ### ridge regression fm <- glmregNB(art ~ ., alpha=0, data = bioChemists, lambda=seq(0.001, 1, by=0.01)) fm <- cv.glmregNB(art ~ ., alpha=0, data = bioChemists, lambda=seq(0.001, 1, by=0.01)) ## End(Not run)## Not run: data("bioChemists", package = "pscl") system.time(fm_nb1 <- glmregNB(art ~ ., data = bioChemists, parallel=FALSE)) system.time(fm_nb2 <- glmregNB(art ~ ., data = bioChemists, parallel=TRUE, n.cores=2)) coef(fm_nb1) ### ridge regression fm <- glmregNB(art ~ ., alpha=0, data = bioChemists, lambda=seq(0.001, 1, by=0.01)) fm <- cv.glmregNB(art ~ ., alpha=0, data = bioChemists, lambda=seq(0.001, 1, by=0.01)) ## End(Not run)
Constructing Hessian matrix for regularized regression parameters.
hessianReg(x, which, ...)hessianReg(x, which, ...)
x |
a fitted model object. |
which |
which penalty parameter(s)? |
... |
arguments passed to the |
hessianReg is a function to compute the Hessian matrix estimate of non-zero regularized estimators. Implemented only for zipath object with family="negbin" in the current version.
A matrix containing the Hessian matrix estimate for the non-zero parameters.
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) hessianReg(fm_zinb, which=which.min(fm_zinb$bic))data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) hessianReg(fm_zinb, which=which.min(fm_zinb$bic))
Fit a robust GLM where the loss function is a composite function cfunodfun.
## S3 method for class 'formula' irglm(formula, data, weights, offset=NULL, contrasts=NULL, cfun="ccave", dfun=gaussian(), s=NULL, delta=0.1, fk=NULL, init.family=NULL, iter=10, reltol=1e-5, theta, x.keep=FALSE, y.keep=TRUE, trace=FALSE, ...)## S3 method for class 'formula' irglm(formula, data, weights, offset=NULL, contrasts=NULL, cfun="ccave", dfun=gaussian(), s=NULL, delta=0.1, fk=NULL, init.family=NULL, iter=10, reltol=1e-5, theta, x.keep=FALSE, y.keep=TRUE, trace=FALSE, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights. |
x |
input matrix, of dimension nobs x nvars; each row is an observation vector |
y |
response variable. Quantitative for |
contrasts |
the contrasts corresponding to |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
cfun |
character, type of convex cap (concave) function.
|
dfun |
character, type of convex component.
|
init.family |
character value for initial family, one of "clossR","closs","gloss","qloss", which can be used to derive an initial estimator, if the selection is different from the default value |
s |
tuning parameter of |
delta |
a small positive number provided by user only if |
fk |
predicted values at an iteration in the IRGLM algorithm |
iter |
number of iteration in the IRGLM algorithm |
reltol |
convergency criteria in the IRGLM algorithm |
theta |
an overdispersion scaling parameter for |
x.keep, y.keep
|
logical values indicating whether the response vector and model matrix used in the fitting process should be returned as components of the returned value, x is a design matrix of dimension n * p, and x is a vector of observations of length n. |
trace |
if |
... |
other arguments passing to |
A robust linear, logistic or Poisson regression model is fit by the iteratively reweighted GLM (IRGLM). The output weights_update is a useful diagnostic to the outlier status of the observations.
An object with S3 class "irglm", "glm" for various types of models.
call |
the call that produced the model fit |
weights |
original weights used in the model |
weights_update |
weights in the final iteration of the IRGLM algorithm |
cfun, s, dfun
|
original input arguments |
is.offset |
is offset used? |
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) fit=irglm(g2~x,data=data.frame(cbind(x, g2)), s=1, cfun="ccave", dfun=gaussian()) fit$weights_updatex=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) fit=irglm(g2~x,data=data.frame(cbind(x, g2)), s=1, cfun="ccave", dfun=gaussian()) fit$weights_update
Fit a robust penalized GLM where the loss function is a composite function cfunodfun + penalty. This is the wrapper function of irglmreg_fit
## S3 method for class 'formula' irglmreg(formula, data, weights, offset=NULL, contrasts=NULL, ...) ## S3 method for class 'matrix' irglmreg(x, y, weights, offset=NULL, ...) ## Default S3 method: irglmreg(x, ...)## S3 method for class 'formula' irglmreg(formula, data, weights, offset=NULL, contrasts=NULL, ...) ## S3 method for class 'matrix' irglmreg(x, y, weights, offset=NULL, ...) ## Default S3 method: irglmreg(x, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights. If |
x |
input matrix, of dimension nobs x nvars; each row is an observation vector |
y |
response variable. Quantitative for |
offset |
Not implemented yet |
contrasts |
the contrasts corresponding to |
... |
Other arguments passing to |
The computing is done by the iteratively reweighted penalized GLM, an application of the iteratively reweighted convex optimization (IRCO). Here convex is the loss function induced by dfun, not the penalty function. The output weights_update is a useful diagnostic to the outlier status of the observations.
The regularization path is computed for the lasso (or elastic net penalty), scad (or snet) and mcp (or mnet penalty), at a grid of values for the regularization parameter lambda.
The sequence of robust models implied by lambda is fit by the IRCO along with coordinate
descent. Note that the objective function is
if standardize=FALSE and
if standardize=TRUE.
An object with S3 class "irglmreg" for the various types of models.
call |
the call that produced this object |
b0 |
Intercept sequence of length |
beta |
A |
lambda |
The actual sequence of |
nobs |
number of observations |
risk |
if |
pll |
if |
fitted.values |
predicted values depending on |
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
print, predict, coef and plot methods, and the cv.irglmreg function.
#binomial x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) fit1=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=TRUE,type.init="bst") #fit1$risk ## Not run: ### different solution paths via a combination of type.path, decreasing and type.init fit1=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=TRUE,type.init="bst") fit2=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=FALSE,type.init="bst") fit3=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="nonactive", decreasing=TRUE,type.init="bst") fit4=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="nonactive", decreasing=FALSE,type.init="bst") fit5=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=TRUE,type.init="co") fit6=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=FALSE,type.init="co") fit7=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="nonactive", decreasing=TRUE,type.init="co") fit8=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="nonactive", decreasing=FALSE,type.init="co") ## End(Not run)#binomial x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) fit1=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=TRUE,type.init="bst") #fit1$risk ## Not run: ### different solution paths via a combination of type.path, decreasing and type.init fit1=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=TRUE,type.init="bst") fit2=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=FALSE,type.init="bst") fit3=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="nonactive", decreasing=TRUE,type.init="bst") fit4=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="nonactive", decreasing=FALSE,type.init="bst") fit5=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=TRUE,type.init="co") fit6=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="active", decreasing=FALSE,type.init="co") fit7=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="nonactive", decreasing=TRUE,type.init="co") fit8=irglmreg(x,g2,s=1,cfun="ccave",dfun="gaussian",type.path="nonactive", decreasing=FALSE,type.init="co") ## End(Not run)
Fit a robust penalized GLM where the loss function is a composite function cfunodfun + penalty. This does computing for irglmreg.
irglmreg_fit(x, y, weights, offset, cfun="ccave", dfun="gaussian", s=NULL, delta=0.1, fk=NULL, iter=10, reltol=1e-5, penalty=c("enet","mnet","snet"), nlambda=100, lambda=NULL, type.path=c("active", "nonactive"), decreasing=TRUE, lambda.min.ratio=ifelse(nobs<nvars,.05, .001), alpha=1, gamma=3, rescale=TRUE, standardize=TRUE, intercept=TRUE, penalty.factor= NULL, maxit=1000, type.init=c("bst", "co", "heu"), init.family=NULL, mstop.init=10, nu.init=0.1, eps=.Machine$double.eps, epscycle=10, thresh=1e-6, parallel=FALSE, n.cores=2, theta, trace=FALSE, tracelevel=1)irglmreg_fit(x, y, weights, offset, cfun="ccave", dfun="gaussian", s=NULL, delta=0.1, fk=NULL, iter=10, reltol=1e-5, penalty=c("enet","mnet","snet"), nlambda=100, lambda=NULL, type.path=c("active", "nonactive"), decreasing=TRUE, lambda.min.ratio=ifelse(nobs<nvars,.05, .001), alpha=1, gamma=3, rescale=TRUE, standardize=TRUE, intercept=TRUE, penalty.factor= NULL, maxit=1000, type.init=c("bst", "co", "heu"), init.family=NULL, mstop.init=10, nu.init=0.1, eps=.Machine$double.eps, epscycle=10, thresh=1e-6, parallel=FALSE, n.cores=2, theta, trace=FALSE, tracelevel=1)
x |
input matrix, of dimension nobs x nvars; each row is an observation vector. |
y |
response variable. Quantitative for |
weights |
observation weights. Can be total counts if responses are proportion matrices. Default is 1 for each observation |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
cfun |
character, type of convex cap (concave) function.
|
dfun |
character, type of convex downward function.
|
s |
tuning parameter of |
delta |
a small positive number provided by user only if |
fk |
predicted values at an iteration in the IRCO algorithm |
nlambda |
The number of |
lambda |
by default, the algorithm provides a sequence of regularization values, or a user supplied |
type.path |
solution path for |
lambda.min.ratio |
Smallest value for |
alpha |
The |
gamma |
The tuning parameter of the |
rescale |
logical value, if TRUE, adaptive rescaling of the penalty parameter for |
standardize |
logical value for x variable standardization, prior to
fitting the model sequence. The coefficients are always returned on
the original scale. Default is |
intercept |
logical value: if TRUE (default), intercept(s) are fitted; otherwise, intercept(s) are set to zero |
penalty.factor |
This is a number that multiplies |
type.init |
a method to determine the initial values. If |
init.family |
character value for initial family, one of "clossR", "closs","gloss","qloss", which can be used to derive an initial estimator, if the selection is different from the default value |
mstop.init |
an integer giving the number of boosting iterations when |
nu.init |
a small number (between 0 and 1) defining the step size or shrinkage parameter when |
decreasing |
only used if |
iter |
number of iteration in the IRCO algorithm |
maxit |
Within each IRCO algorithm iteration, maximum number of coordinate descent iterations for each |
reltol |
convergency criteria in the IRCO algorithm |
eps |
If a coefficient is less than |
epscycle |
If |
thresh |
Convergence threshold for coordinate descent. Defaults value is |
penalty |
Type of regularization |
theta |
an overdispersion scaling parameter for |
parallel, n.cores
|
If |
trace, tracelevel
|
If |
A case weighted penalized least squares or GLM is fit by the iteratively reweighted convex optimization (IRCO), where the loss function is a composite function cfunodfun + penalty. Here convex is the loss function induced by dfun, not the penalty function.
The sequence of robust models implied by lambda is fit by IRCO along with coordinate
descent. Note that the objective function is
if standardize=FALSE and
if standardize=TRUE.
An object with S3 class "irglmreg" for the various types of models.
call |
the call that produced the model fit |
b0 |
Intercept sequence of length |
beta |
A |
lambda |
The actual sequence of |
weights_update |
A |
decreasing |
if |
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
Fit case weighted support vector machines with robust loss functions. This is the wrapper function of irsvm_fit, which does the computing.
## S3 method for class 'formula' irsvm(formula, data, weights, contrasts=NULL, ...) ## S3 method for class 'matrix' irsvm(x, y, weights, ...) ## Default S3 method: irsvm(x, ...)## S3 method for class 'formula' irsvm(formula, data, weights, contrasts=NULL, ...) ## S3 method for class 'matrix' irsvm(x, y, weights, ...) ## Default S3 method: irsvm(x, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights |
x |
input matrix, of dimension nobs x nvars; each row is an observation vector |
y |
response variable. Quantitative for |
contrasts |
the contrasts corresponding to |
... |
Other arguments passing to |
Fit a robust SVM where the loss function is a composite function cfunotype + penalty.
The model is fit by the iteratively reweighted SVM, an application of the iteratively reweighted convex optimization (IRCO). Here convex is the loss function induced by type.
For linear kernel, the coefficients of the regression/decision hyperplane
can be extracted using the coef method.
An object with S3 class "wsvm" for various types of models.
call |
the call that produced this object |
weights_update |
weights in the final iteration of the IRCO algorithm |
cfun, s
|
original input arguments |
delta |
delta value used for |
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
irsvm_fit, print, predict, coef.
#binomial x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) fit=irsvm(x,g2,s=1,cfun="ccave",type="C-classification")#binomial x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) fit=irsvm(x,g2,s=1,cfun="ccave",type="C-classification")
irsvm_fit is used to train a subject weighted support vector machine where the weights are provided iteratively from robust loss function with the iteratively reweighted convex optimization (IRCO). It can be used to carry out robust regression and binary classification. This does computing for the wrapper function irsvm.
irsvm_fit(x, y, weights, cfun="ccave", s=NULL, delta=0.0001, type = NULL, kernel="radial", cost=1, epsilon = 0.1, iter=10, reltol=1e-5, trace=FALSE, ...)irsvm_fit(x, y, weights, cfun="ccave", s=NULL, delta=0.0001, type = NULL, kernel="radial", cost=1, epsilon = 0.1, iter=10, reltol=1e-5, trace=FALSE, ...)
x |
a data matrix, a vector, or a sparse 'design matrix' (object of class
|
y |
a response vector with one label for each row/component of
|
weights |
the weight of each subject. It should be in the same length of |
cfun |
character, type of convex cap (concave) function.
|
s |
tuning parameter of |
delta |
a small positive number provided by user only if |
type |
|
kernel |
the kernel used in training and predicting. You
might consider changing some of the following parameters, depending
on the kernel type.
|
cost |
cost of constraints violation (default: 1)—it is the
‘C’-constant of the regularization term in the Lagrange formulation. This is proportional to the inverse of |
epsilon |
epsilon in the insensitive-loss function (default: 0.1) |
iter |
number of iteration in the IRCO algorithm |
reltol |
convergency criteria in the IRCO algorithm |
trace |
If |
... |
additional parameters for function |
A case weighted SVM is fit by the IRCO algorithm, where the loss function is a composite function of cfunotype, plus a penalty.
Additional arguments include degree, gamma, coef0,
class.weights, cachesize, tolerance, shrinking, propbability, fitted, the same as "wsvm" in package WeightSVM.
An object of class "wsvm" (see package WeightSVM) containing the fitted model, including:
SV |
The resulting support vectors (possibly scaled). |
index |
The index of the resulting support vectors in the data
matrix. Note that this index refers to the preprocessed data (after
the possible effect of |
coefs |
The corresponding coefficients times the training labels. |
rho |
The negative intercept. |
sigma |
In case of a probabilistic regression model, the scale parameter of the hypothesized (zero-mean) laplace distribution estimated by maximum likelihood. |
probA, probB
|
numeric vectors of length 2, number of classes, containing the parameters of the logistic distributions fitted to the decision values of the binary classifiers (1 / (1 + exp(a x + b))). |
Zhu Wang [email protected]
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
irsvm, print, predict, coef and plot.
data(iris) iris <- subset(iris, Species %in% c("setosa", "versicolor")) # default with factor response: model <- irsvm(Species ~ ., data = iris, kernel="linear", trace=TRUE) model <- irsvm(Species ~ ., data = iris) # alternatively the traditional interface: x <- subset(iris, select = -Species) y <- iris$Species model <- irsvm(x, y) # test with train data pred <- predict(model, x) # (same as:) pred <- fitted(model) # Check accuracy: table(pred, y) # compute decision values and probabilities: pred <- predict(model, x, decision.values = TRUE) attr(pred, "decision.values") # visualize (classes by color, SV by crosses): plot(cmdscale(dist(iris[,-5])), col = as.integer(iris[,5]), pch = c("o","+")[1:100 %in% model$index + 1]) ## try regression mode on two dimensions # create data x <- seq(0.1, 5, by = 0.05) y <- log(x) + rnorm(x, sd = 0.2) # estimate model and predict input values m <- irsvm(x, y) new <- predict(m, x) # visualize plot(x, y) points(x, log(x), col = 2) points(x, new, col = 4)data(iris) iris <- subset(iris, Species %in% c("setosa", "versicolor")) # default with factor response: model <- irsvm(Species ~ ., data = iris, kernel="linear", trace=TRUE) model <- irsvm(Species ~ ., data = iris) # alternatively the traditional interface: x <- subset(iris, select = -Species) y <- iris$Species model <- irsvm(x, y) # test with train data pred <- predict(model, x) # (same as:) pred <- fitted(model) # Check accuracy: table(pred, y) # compute decision values and probabilities: pred <- predict(model, x, decision.values = TRUE) attr(pred, "decision.values") # visualize (classes by color, SV by crosses): plot(cmdscale(dist(iris[,-5])), col = as.integer(iris[,5]), pch = c("o","+")[1:100 %in% model$index + 1]) ## try regression mode on two dimensions # create data x <- seq(0.1, 5, by = 0.05) y <- log(x) + rnorm(x, sd = 0.2) # estimate model and predict input values m <- irsvm(x, y) new <- predict(m, x) # visualize plot(x, y) points(x, log(x), col = 2) points(x, new, col = 4)
Compute composite loss value
loss2(y, f, weights, cfun, dfun, s, delta=0.0001)loss2(y, f, weights, cfun, dfun, s, delta=0.0001)
y |
response variable values |
f |
linear predictor values of |
weights |
observation weights, same length as |
cfun |
integer from 1-8, concave function as in |
dfun |
integer from 1-7, convex function as in |
s |
tuning parameter of |
delta |
a small positive number provided by user only if |
An internal function. For large s values, the loss can be 0 with cfun=2,3,4, or "acave", "bcave", "ccave".
Weighted loss values
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
loss3
irglm
irglmreg
loss2_irsvm
Compute composite loss value for epsilon-insensitive type function
loss2_irsvm(y, f, weights, cfun, dfun, s, eps, delta=0.0001)loss2_irsvm(y, f, weights, cfun, dfun, s, eps, delta=0.0001)
y |
response variable values |
f |
fitted values of |
weights |
observation weights, same length as |
cfun |
integer from 1-8, concave function as in |
dfun |
integer value, only |
s |
tuning parameter of |
delta |
a small positive number provided by user only if |
eps |
non-negative parameter for epsilon-insensitive loss |
For large s values, the loss can be 0 with cfun=2,3,4, or "acave", "bcave", "ccave".
Weighted loss values
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
Compute composite loss value
loss3(y, mu, theta, weights, cfun, family, s, delta)loss3(y, mu, theta, weights, cfun, family, s, delta)
y |
response variable values, 0/1 if |
mu |
response prediction of |
theta |
scale parameter for |
weights |
observation weights, same length as |
cfun |
integer from 1-8, concave function as in |
family |
integer 2, 3 or 4, convex function binomial, Poisson or negative binomial, respectively |
s |
tuning parameter of |
delta |
a small positive number provided by user only if |
For large s values, the loss can be 0 with cfun=2,3,4, or "acave", "bcave", "ccave".
Weighted loss values
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
loss2
irglm
irglmreg
loss2_irsvm
Estimating the variance of the first derivative of log-likelihood function
meatReg(x, which, ...)meatReg(x, which, ...)
x |
a fitted model object. Currently only implemented for |
which |
which penalty parameter(s)? |
... |
arguments passed to the |
See reference below
A
covariance matrix of first derivative of log-likelihood function
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
sandwichReg, breadReg, estfunReg
data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) meatReg(fm_zinb, which=which.min(fm_zinb$bic))data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) meatReg(fm_zinb, which=which.min(fm_zinb$bic))
Methods for models fitted by coordinate descent algorithms.
## S3 method for class 'glmreg' AIC(object, ..., k) ## S3 method for class 'zipath' AIC(object, ..., k) ## S3 method for class 'glmreg' BIC(object, ...) ## S3 method for class 'zipath' BIC(object, ...)## S3 method for class 'glmreg' AIC(object, ..., k) ## S3 method for class 'zipath' AIC(object, ..., k) ## S3 method for class 'glmreg' BIC(object, ...) ## S3 method for class 'zipath' BIC(object, ...)
object |
objects of class |
... |
additional arguments passed to calls. |
k |
numeric, the penalty per parameter to be used; the default
|
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
Fit a linear model via penalized nonconvex loss function.
## S3 method for class 'formula' ncl(formula, data, weights, offset=NULL, contrasts=NULL, x.keep=FALSE, y.keep=TRUE, ...) ## S3 method for class 'matrix' ncl(x, y, weights, offset=NULL, ...) ## Default S3 method: ncl(x, ...)## S3 method for class 'formula' ncl(formula, data, weights, offset=NULL, contrasts=NULL, x.keep=FALSE, y.keep=TRUE, ...) ## S3 method for class 'matrix' ncl(x, y, weights, offset=NULL, ...) ## Default S3 method: ncl(x, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights. If |
x |
input matrix, of dimension nobs x nvars; each row is an observation vector |
y |
response variable. Quantitative for |
offset |
Not implemented yet |
contrasts |
the contrasts corresponding to |
x.keep, y.keep
|
For glmreg: logical values indicating whether the response vector and model matrix used in the fitting process should be returned as components of the returned value. For ncl_fit: x is a design matrix of dimension n * p, and x is a vector of observations of length n. |
... |
Other arguments passing to |
The robust linear model is fit by majorization-minimization along with linear regression. Note that the objective function is
.
An object with S3 class "ncl" for the various types of models.
call |
the call that produced this object |
fitted.values |
predicted values |
h |
pseudo response values in the MM algorithm |
Zhu Wang <[email protected]>
Zhu Wang (2021), MM for Penalized Estimation, TEST, doi:10.1007/s11749-021-00770-2
#binomial x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) fit=ncl(x,g2,s=1,rfamily="closs")#binomial x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) fit=ncl(x,g2,s=1,rfamily="closs")
Fit a linear model via penalized nonconvex loss function.
ncl_fit(x,y, weights, offset=NULL, rfamily=c("clossR", "closs", "gloss", "qloss"), s=NULL, fk=NULL, iter=10, reltol=1e-5, trace=FALSE)ncl_fit(x,y, weights, offset=NULL, rfamily=c("clossR", "closs", "gloss", "qloss"), s=NULL, fk=NULL, iter=10, reltol=1e-5, trace=FALSE)
x |
input matrix, of dimension nobs x nvars; each row is an observation vector. |
y |
response variable. Quantitative for |
weights |
observation weights. Can be total counts if responses are proportion matrices. Default is 1 for each observation |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
rfamily |
Response type and relevant loss functions (see above) |
s |
nonconvex loss tuning parameter for robust regression and classification. |
fk |
predicted values at an iteration in the MM algorithm |
iter |
number of iteration in the MM algorithm |
reltol |
convergency criteria |
trace |
If |
The robust linear model is fit by majorization-minimization along with least squares. Note that the objective function is
.
An object with S3 class "ncl" for the various types of models.
call |
the call that produced the model fit |
fitted.values |
predicted values |
h |
pseudo response values in the MM algorithm |
Zhu Wang <[email protected]>
Zhu Wang (2021), MM for Penalized Estimation, TEST, doi:10.1007/s11749-021-00770-2
Fit a linear model via penalized nonconvex loss function. The regularization path is computed for the lasso (or elastic net penalty), scad (or snet) and mcp (or mnet penalty), at a grid of values for the regularization parameter lambda. The name refers to NonConvex Loss with REGularization.
## S3 method for class 'formula' nclreg(formula, data, weights, offset=NULL, contrasts=NULL, ...) ## S3 method for class 'matrix' nclreg(x, y, weights, offset=NULL, ...) ## Default S3 method: nclreg(x, ...)## S3 method for class 'formula' nclreg(formula, data, weights, offset=NULL, contrasts=NULL, ...) ## S3 method for class 'matrix' nclreg(x, y, weights, offset=NULL, ...) ## Default S3 method: nclreg(x, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights. If |
x |
input matrix, of dimension nobs x nvars; each row is an observation vector |
y |
response variable. Quantitative for |
offset |
Not implemented yet |
contrasts |
the contrasts corresponding to |
... |
Other arguments passing to |
The sequence of robust models implied by lambda is fit by majorization-minimization along with coordinate
descent. Note that the objective function is
if standardize=FALSE and
if standardize=TRUE.
An object with S3 class "nclreg" for the various types of models.
call |
the call that produced this object |
b0 |
Intercept sequence of length |
beta |
A |
lambda |
The actual sequence of |
nobs |
number of observations |
risk |
if |
pll |
if |
fitted.values |
predicted values depending on |
Zhu Wang <[email protected]>
Zhu Wang (2021), MM for Penalized Estimation, TEST, doi:10.1007/s11749-021-00770-2
print, predict, coef and plot methods, and the cv.nclreg function.
#binomial x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) ### different solution paths via a combination of type.path, decreasing and type.init fit1=nclreg(x,g2,s=1,rfamily="closs",type.path="active",decreasing=TRUE,type.init="bst") fit2=nclreg(x,g2,s=1,rfamily="closs",type.path="active",decreasing=FALSE,type.init="bst") fit3=nclreg(x,g2,s=1,rfamily="closs",type.path="nonactive",decreasing=TRUE,type.init="bst") fit4=nclreg(x,g2,s=1,rfamily="closs",type.path="nonactive",decreasing=FALSE,type.init="bst") fit5=nclreg(x,g2,s=1,rfamily="closs",type.path="active",decreasing=TRUE,type.init="ncl") fit6=nclreg(x,g2,s=1,rfamily="closs",type.path="active",decreasing=FALSE,type.init="ncl") fit7=nclreg(x,g2,s=1,rfamily="closs",type.path="nonactive",decreasing=TRUE,type.init="ncl") fit8=nclreg(x,g2,s=1,rfamily="closs",type.path="nonactive",decreasing=FALSE,type.init="ncl")#binomial x=matrix(rnorm(100*20),100,20) g2=sample(c(-1,1),100,replace=TRUE) ### different solution paths via a combination of type.path, decreasing and type.init fit1=nclreg(x,g2,s=1,rfamily="closs",type.path="active",decreasing=TRUE,type.init="bst") fit2=nclreg(x,g2,s=1,rfamily="closs",type.path="active",decreasing=FALSE,type.init="bst") fit3=nclreg(x,g2,s=1,rfamily="closs",type.path="nonactive",decreasing=TRUE,type.init="bst") fit4=nclreg(x,g2,s=1,rfamily="closs",type.path="nonactive",decreasing=FALSE,type.init="bst") fit5=nclreg(x,g2,s=1,rfamily="closs",type.path="active",decreasing=TRUE,type.init="ncl") fit6=nclreg(x,g2,s=1,rfamily="closs",type.path="active",decreasing=FALSE,type.init="ncl") fit7=nclreg(x,g2,s=1,rfamily="closs",type.path="nonactive",decreasing=TRUE,type.init="ncl") fit8=nclreg(x,g2,s=1,rfamily="closs",type.path="nonactive",decreasing=FALSE,type.init="ncl")
Fit a linear model via penalized nonconvex loss function. The regularization path is computed for the lasso (or elastic net penalty), scad (or snet) and mcp (or mnet penalty), at a grid of values for the regularization parameter lambda.
nclreg_fit(x, y, weights, offset, rfamily=c("clossR", "closs", "gloss", "qloss"), s=NULL, fk=NULL, iter=10, reltol=1e-5, penalty=c("enet","mnet","snet"), nlambda=100,lambda=NULL, type.path=c("active", "nonactive", "onestep"), decreasing=FALSE, lambda.min.ratio=ifelse(nobs<nvars,.05, .001), alpha=1, gamma=3, standardize=TRUE, intercept=TRUE, penalty.factor=NULL, maxit=1000, type.init=c("bst", "ncl", "heu"), mstop.init=10, nu.init=0.1, eps=.Machine$double.eps, epscycle=10, thresh=1e-6, trace=FALSE)nclreg_fit(x, y, weights, offset, rfamily=c("clossR", "closs", "gloss", "qloss"), s=NULL, fk=NULL, iter=10, reltol=1e-5, penalty=c("enet","mnet","snet"), nlambda=100,lambda=NULL, type.path=c("active", "nonactive", "onestep"), decreasing=FALSE, lambda.min.ratio=ifelse(nobs<nvars,.05, .001), alpha=1, gamma=3, standardize=TRUE, intercept=TRUE, penalty.factor=NULL, maxit=1000, type.init=c("bst", "ncl", "heu"), mstop.init=10, nu.init=0.1, eps=.Machine$double.eps, epscycle=10, thresh=1e-6, trace=FALSE)
x |
input matrix, of dimension nobs x nvars; each row is an observation vector. |
y |
response variable. Quantitative for |
weights |
observation weights. Can be total counts if responses are proportion matrices. Default is 1 for each observation |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. Currently only one offset term can be included in the formula. |
rfamily |
Response type and relevant loss functions (see above) |
s |
nonconvex loss tuning parameter for robust regression and classification. The |
fk |
predicted values at an iteration in the MM algorithm |
nlambda |
The number of |
lambda |
by default, the algorithm provides a sequence of regularization values, or a user supplied |
type.path |
solution path. If |
lambda.min.ratio |
Smallest value for |
alpha |
The |
gamma |
The tuning parameter of the |
standardize |
logical value for x variable standardization, prior to
fitting the model sequence. The coefficients are always returned on
the original scale. Default is |
intercept |
logical value: if TRUE (default), intercept(s) are fitted; otherwise, intercept(s) are set to zero |
penalty.factor |
This is a number that multiplies |
type.init |
a method to determine the initial values. If |
mstop.init |
an integer giving the number of boosting iterations when |
nu.init |
a small number (between 0 and 1) defining the step size or shrinkage parameter when |
decreasing |
only used if |
iter |
number of iteration in the MM algorithm |
maxit |
Within each MM algorithm iteration, maximum number of coordinate descent iterations for each |
reltol |
convergency criteria |
eps |
If a coefficient is less than |
epscycle |
If |
thresh |
Convergence threshold for coordinate descent. Defaults value is |
penalty |
Type of regularization |
trace |
If |
The sequence of robust models implied by lambda is fit by majorization-minimization along with coordinate
descent. Note that the objective function is
if standardize=FALSE and
if standardize=TRUE.
An object with S3 class "nclreg" for the various types of models.
call |
the call that produced the model fit |
b0 |
Intercept sequence of length |
beta |
A |
lambda |
The actual sequence of |
decreasing |
if |
Zhu Wang <[email protected]>
Zhu Wang (2021), MM for Penalized Estimation, TEST, doi:10.1007/s11749-021-00770-2
Produces a coefficient profile plot of the coefficient paths for a
fitted "glmreg" object.
## S3 method for class 'glmreg' plot(x, xvar = c("norm", "lambda", "dev"), label = FALSE, shade=TRUE, ...)## S3 method for class 'glmreg' plot(x, xvar = c("norm", "lambda", "dev"), label = FALSE, shade=TRUE, ...)
x |
fitted |
xvar |
What is on the X-axis. |
label |
If |
shade |
Should nonconvex region be shaded? Default is TRUE. Code developed for all |
... |
Other graphical parameters to plot |
A coefficient profile plot is produced.
Zhu Wang [email protected]
glmreg, and print, predict and coef methods.
x=matrix(rnorm(100*20),100,20) y=rnorm(100) fit1=glmreg(x,y) plot(fit1) plot(fit1,xvar="lambda",label=TRUE)x=matrix(rnorm(100*20),100,20) y=rnorm(100) fit1=glmreg(x,y) plot(fit1) plot(fit1,xvar="lambda",label=TRUE)
This function returns predictions from
a fitted "glmreg" object.
## S3 method for class 'glmreg' predict(object,newx,newoffset,which=1:length(object$lambda), type=c("link","response","class","coefficients","nonzero"), na.action=na.pass, ...) ## S3 method for class 'glmreg' coef(object,which=1:length(object$lambda),...)## S3 method for class 'glmreg' predict(object,newx,newoffset,which=1:length(object$lambda), type=c("link","response","class","coefficients","nonzero"), na.action=na.pass, ...) ## S3 method for class 'glmreg' coef(object,which=1:length(object$lambda),...)
object |
Fitted |
newx |
Matrix of values at which predictions are to be made. Not
used for |
which |
Indices of the penalty parameter |
type |
Type of prediction: |
newoffset |
an offset term used in prediction |
na.action |
action for missing data value |
... |
arguments for predict |
The returned object depends on type.
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
## Dobson (1990) Page 93: Randomized Controlled Trial : counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9) treatment <- gl(3,3) print(d.AD <- data.frame(treatment, outcome, counts)) fit <- glmreg(counts ~ outcome + treatment, data=d.AD, family="poisson") predict(fit, newx=d.AD[,1:2]) summary(fit) coef(fit)## Dobson (1990) Page 93: Randomized Controlled Trial : counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9) treatment <- gl(3,3) print(d.AD <- data.frame(treatment, outcome, counts)) fit <- glmreg(counts ~ outcome + treatment, data=d.AD, family="poisson") predict(fit, newx=d.AD[,1:2]) summary(fit) coef(fit)
Methods for extracting information from fitted penalized zero-inflated
regression model objects of class "zipath".
## S3 method for class 'zipath' predict(object, newdata, which = 1:object$nlambda, type = c("response", "prob", "count", "zero", "nonzero"), na.action = na.pass, at = NULL, ...) ## S3 method for class 'zipath' residuals(object, type = c("pearson", "response"), ...) ## S3 method for class 'zipath' coef(object, which=1:object$nlambda, model = c("full", "count", "zero"), ...) ## S3 method for class 'zipath' terms(x, model = c("count", "zero"), ...) ## S3 method for class 'zipath' model.matrix(object, model = c("count", "zero"), ...)## S3 method for class 'zipath' predict(object, newdata, which = 1:object$nlambda, type = c("response", "prob", "count", "zero", "nonzero"), na.action = na.pass, at = NULL, ...) ## S3 method for class 'zipath' residuals(object, type = c("pearson", "response"), ...) ## S3 method for class 'zipath' coef(object, which=1:object$nlambda, model = c("full", "count", "zero"), ...) ## S3 method for class 'zipath' terms(x, model = c("count", "zero"), ...) ## S3 method for class 'zipath' model.matrix(object, model = c("count", "zero"), ...)
object, x
|
an object of class |
newdata |
optionally, a data frame in which to look for variables with which to predict. If omitted, the original observations are used. |
which |
Indices of the penalty parameters |
type |
character specifying the type of predictions or residuals, respectively. For details see below. |
na.action |
function determining what should be done with missing values
in |
at |
optionally, if |
model |
character specifying for which component of the model the terms or model matrix should be extracted. |
... |
currently not used. |
Re-uses the design of function zeroinfl in package pscl (see reference). A set of standard extractor functions for fitted model objects is available for
objects of class "zipath", including methods to the generic functions
print and summary which print the estimated
coefficients along with some further information.
As usual, the summary method returns an object of class "summary.zipath"
containing the relevant summary statistics which can subsequently be printed
using the associated print method.
The methods for coef by default
return a single vector of coefficients and their associated covariance matrix,
respectively, i.e., all coefficients are concatenated. By setting the model
argument, the estimates for the corresponding model components can be extracted.
Both the fitted and predict methods can
compute fitted responses. The latter additionally provides the predicted density
(i.e., probabilities for the observed counts), the predicted mean from the count
component (without zero inflation) and the predicted probability for the zero
component. The residuals method can compute
raw residuals (observed - fitted) and Pearson residuals (raw residuals scaled by
square root of variance function).
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
## Not run: data("bioChemists", package = "pscl") fm_zip <- zipath(art ~ . | ., data = bioChemists, nlambda=10) plot(residuals(fm_zip) ~ fitted(fm_zip)) coef(fm_zip, model = "count") coef(fm_zip, model = "zero") summary(fm_zip) logLik(fm_zip) ## End(Not run)## Not run: data("bioChemists", package = "pscl") fm_zip <- zipath(art ~ . | ., data = bioChemists, nlambda=10) plot(residuals(fm_zip) ~ fitted(fm_zip)) coef(fm_zip, model = "count") coef(fm_zip, model = "zero") summary(fm_zip) logLik(fm_zip) ## End(Not run)
compute p-values from penalized zero-inflated Poisson, negative binomial and geometric model with multi-split data
pval.zipath(formula, data, weights, subset, na.action, offset, standardize=TRUE, family = c("poisson", "negbin", "geometric"), penalty = c("enet", "mnet", "snet"), gamma.count = 3, gamma.zero = 3, prop=0.5, trace=TRUE, B=10, ...)pval.zipath(formula, data, weights, subset, na.action, offset, standardize=TRUE, family = c("poisson", "negbin", "geometric"), penalty = c("enet", "mnet", "snet"), gamma.count = 3, gamma.zero = 3, prop=0.5, trace=TRUE, B=10, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights. If |
subset |
subset of data |
na.action |
how to deal with missing data |
offset |
Not implemented yet |
standardize |
logical value, should variables be standardized? |
family |
family to fit |
penalty |
penalty considered as one of |
gamma.count |
The tuning parameter of the |
gamma.zero |
The tuning parameter of the |
prop |
proportion of data split, default is 50/50 split |
trace |
logical value, if TRUE, print detailed calculation results |
B |
number of repeated multi-split replications |
... |
Other arguments passing to |
compute p-values from penalized zero-inflated Poisson, negative binomial and geometric model with multi-split data
count.pval |
raw p-values in the count component |
zero.pval |
raw p-values in the zero component |
count.pval.q |
Q value for the count component |
zero.pval.q |
Q value for the zero component |
Zhu Wang <[email protected]>
Nicolai Meinshausen, Lukas Meier and Peter Buehlmann (2013) p-Values for High-Dimensional Regression, Journal of the American Statistical Association, 104(488), 1671–1681.
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
random number generation of zero-inflated count response
rzi(n, x, z, a, b, theta=1, family=c("poisson", "negbin", "geometric"), infl=TRUE)rzi(n, x, z, a, b, theta=1, family=c("poisson", "negbin", "geometric"), infl=TRUE)
n |
sample size of random number generation |
x |
design matrix of count model |
z |
design matrix of zero model |
a |
coefficient vector for |
b |
coefficient vector for |
theta |
dispersion parameter for |
family |
distribution of count model |
infl |
logical value, if TRUE, zero-inflated count response |
random number generation of zero-inflated count response
numeric vector of zero-inflated count response
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
Constructing sandwich covariance matrix estimators by multiplying bread and meat matrices for regularized regression parameters.
sandwichReg(x, breadreg.=breadReg, meatreg.=meatReg, which, log=FALSE, ...)sandwichReg(x, breadreg.=breadReg, meatreg.=meatReg, which, log=FALSE, ...)
x |
a fitted model object. |
breadreg. |
either a breadReg matrix or a function for computing
this via |
meatreg. |
either a breadReg matrix or a function for computing
this via |
which |
which penalty parameters(s) to compute? |
log |
if TRUE, the corresponding element is with respect to log(theta) in negative binomial regression. Otherwise, for theta |
... |
arguments passed to the |
sandwichReg is a function to compute an estimator for the covariance of the non-zero parameters. It takes a breadReg matrix (i.e., estimator of the expectation of the negative
derivative of the penalized estimating functions) and a meatReg matrix (i.e.,
estimator of the variance of the log-likelihood function) and multiplies
them to a sandwich with meat between two slices of bread. By default
breadReg and meatReg are called. Implemented only for zipath object with family="negbin" in the current version.
A matrix containing the sandwich covariance matrix estimate for the non-zero parameters.
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) sandwichReg(fm_zinb, which=which.min(fm_zinb$bic))data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) sandwichReg(fm_zinb, which=which.min(fm_zinb$bic))
Generic function for computing standard errors of non-zero regularized estimators
se(x, which, log=TRUE, ...)se(x, which, log=TRUE, ...)
x |
a fitted model object. |
which |
which penalty parameter(s)? |
log |
if TRUE, the computed standard error is for log(theta) for negative binomial regression, otherwise, for theta. |
... |
arguments passed to methods. |
A vector containing standard errors of non-zero regularized estimators.
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) res <- se(fm_zinb, which=which.min(fm_zinb$bic))data("bioChemists", package = "pscl") fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10, maxit.em=1) res <- se(fm_zinb, which=which.min(fm_zinb$bic))
Standardize variables. For each column, return mean 0 and mean value of sum of squares = 1.
stan(x, weights)stan(x, weights)
x |
numeric variables, can be a matrix or vector |
weights |
numeric positive vector of weights |
A list with the following items.
x |
standardized variables with each column: mean value 0 and mean value of sum of squares = 1. |
meanx |
a vector of means for each column in the original |
normx |
a vector of scales for each column in the original |
Zhu Wang <[email protected]>
Summary results of fitted penalized negative binomial regression model
## S3 method for class 'glmregNB' summary(object, ...)## S3 method for class 'glmregNB' summary(object, ...)
object |
fitted model object of class |
... |
arguments passed to or from other methods. |
This function is a method for the generic function
summary() for class "glmregNB".
It can be invoked by calling summary(x) for an
object x of the appropriate class, or directly by
calling summary.glmregNB(x) regardless of the
class of the object.
Summary of fitted penalized negative binomial model
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Fit penalized zero-inflated models, generate multiple paths with varying penalty parameters, therefore determine optimal path with respect to a particular penalty parameter
tuning.zipath(formula, data, weights, subset, na.action, offset, standardize=TRUE, family = c("poisson", "negbin", "geometric"), penalty = c("enet", "mnet", "snet"), lambdaCountRatio = .0001, lambdaZeroRatio = c(.1, .01, .001), maxit.theta=1, gamma.count=3, gamma.zero=3, ...)tuning.zipath(formula, data, weights, subset, na.action, offset, standardize=TRUE, family = c("poisson", "negbin", "geometric"), penalty = c("enet", "mnet", "snet"), lambdaCountRatio = .0001, lambdaZeroRatio = c(.1, .01, .001), maxit.theta=1, gamma.count=3, gamma.zero=3, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights. If |
subset |
subset of data |
na.action |
how to deal with missing data |
offset |
Not implemented yet |
standardize |
logical value, should variables be standardized? |
family |
family to fit |
penalty |
penalty considered as one of |
lambdaCountRatio, lambdaZeroRatio
|
Smallest value for |
maxit.theta |
For family="negbin", the maximum iteration allowed for estimating scale parameter theta. Note, the default value 1 is for computing speed purposes, and is typically too small and less desirable in real data analysis |
gamma.count |
The tuning parameter of the |
gamma.zero |
The tuning parameter of the |
... |
Other arguments passing to |
From the default lambdaZeroRatio = c(.1, .01, .001) values,
find optimal lambdaZeroRatio for penalized zero-inflated Poisson, negative binomial and geometric model
An object of class zipath with the optimal lambdaZeroRatio
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
## Not run: ## data data("bioChemists", package = "pscl") ## inflation with regressors ## ("art ~ . | ." is "art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment") fm_zip2 <- tuning.zipath(art ~ . | ., data = bioChemists, nlambda=10) summary(fm_zip2) fm_zinb2 <- tuning.zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) summary(fm_zinb2) ## End(Not run)## Not run: ## data data("bioChemists", package = "pscl") ## inflation with regressors ## ("art ~ . | ." is "art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment") fm_zip2 <- tuning.zipath(art ~ . | ., data = bioChemists, nlambda=10) summary(fm_zip2) fm_zinb2 <- tuning.zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) summary(fm_zinb2) ## End(Not run)
Compute weight value
update_wt(y, ypre, weights, cfun, s, dfun, delta=0.0001)update_wt(y, ypre, weights, cfun, s, dfun, delta=0.0001)
y |
input value of response variable |
ypre |
predicted value of response variable |
weights |
optional numeric vector of weights. |
cfun |
integer from 1-8, concave function as in |
dfun |
integer value, convex function as in |
s |
a numeric value, see details in |
delta |
a positive small value, see details in |
Weight value
Zhu Wang <[email protected]>
Zhu Wang (2024) Unified Robust Estimation, Australian & New Zealand Journal of Statistics. 66(1):77-102.
Fit zero-inflated regression models for count data via penalized maximum likelihood.
## S3 method for class 'formula' zipath(formula, data, weights, offset=NULL, contrasts=NULL, ... ) ## S3 method for class 'matrix' zipath(X, Z, Y, weights, offsetx=NULL, offsetz=NULL, ...) ## Default S3 method: zipath(X, ...)## S3 method for class 'formula' zipath(formula, data, weights, offset=NULL, contrasts=NULL, ... ) ## S3 method for class 'matrix' zipath(X, Z, Y, weights, offsetx=NULL, offsetz=NULL, ...) ## Default S3 method: zipath(X, ...)
formula |
symbolic description of the model, see details. |
data |
argument controlling formula processing
via |
weights |
optional numeric vector of weights. |
offset |
optional numeric vector with an a priori known component to be included in the linear predictor of the count model or zero model. See below for an example. |
contrasts |
a list with elements |
X |
predictor matrix of the count model |
Z |
predictor matrix of the zero model |
Y |
response variable |
offsetx, offsetz
|
optional numeric vector with an a priori known component to be included in the linear predictor of the count model (offsetx)or zero model (offsetz). |
... |
Other arguments which can be passed to |
An object of class "zipath", i.e., a list with components including
coefficients |
a list with elements |
residuals |
a vector of raw residuals (observed - fitted), |
fitted.values |
a vector of fitted means, |
weights |
the case weights used, |
terms |
a list with elements |
theta |
estimate of the additional |
loglik |
log-likelihood of the fitted model, |
family |
character string describing the count distribution used, |
link |
character string describing the link of the zero-inflation model, |
linkinv |
the inverse link function corresponding to |
converged |
logical value, TRUE indicating successful convergence of |
call |
the original function call |
formula |
the original formula |
levels |
levels of the categorical regressors |
contrasts |
a list with elements |
model |
the full model frame (if |
y |
the response count vector (if |
x |
a list with elements |
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.
## data data("bioChemists", package = "pscl") ## with simple inflation (no regressors for zero component) fm_zip <- zipath(art ~ 1 | ., data = bioChemists, nlambda=10) summary(fm_zip) fm_zip <- zipath(art ~ . | 1, data = bioChemists, nlambda=10) summary(fm_zip) ## Not run: fm_zip <- zipath(art ~ . | 1, data = bioChemists, nlambda=10) summary(fm_zip) fm_zinb <- zipath(art ~ . | 1, data = bioChemists, family = "negbin", nlambda=10) summary(fm_zinb) ## inflation with regressors ## ("art ~ . | ." is "art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment") fm_zip2 <- zipath(art ~ . | ., data = bioChemists, nlambda=10) summary(fm_zip2) fm_zinb2 <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) summary(fm_zinb2) ### non-penalized regression, compare with zeroinfl fm_zinb3 <- zipath(art ~ . | ., data = bioChemists, family = "negbin", lambda.count=0, lambda.zero=0, reltol=1e-12) summary(fm_zinb3) library("pscl") fm_zinb4 <- zeroinfl(art ~ . | ., data = bioChemists, dist = "negbin") summary(fm_zinb4) ### offset exposure <- rep(0.5, dim(bioChemists)[1]) fm_zinb <- zipath(art ~ . +offset(log(exposure))| ., data = bioChemists, family = "poisson", nlambda=10) coef <- coef(fm_zinb) ### offset can't be specified in predict function as it has been contained pred <- predict(fm_zinb) ## without inflation ## ("art ~ ." is "art ~ fem + mar + kid5 + phd + ment") fm_pois <- glmreg(art ~ ., data = bioChemists, family = "poisson") coef <- coef(fm_pois) fm_nb <- glmregNB(art ~ ., data = bioChemists) coef <- coef(fm_nb) ### high-dimensional #R CMD check --use-valgrind can be too time extensive for the following model #bioChemists <- cbind(matrix(rnorm(915*100), nrow=915), bioChemists) #fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) ## End(Not run)## data data("bioChemists", package = "pscl") ## with simple inflation (no regressors for zero component) fm_zip <- zipath(art ~ 1 | ., data = bioChemists, nlambda=10) summary(fm_zip) fm_zip <- zipath(art ~ . | 1, data = bioChemists, nlambda=10) summary(fm_zip) ## Not run: fm_zip <- zipath(art ~ . | 1, data = bioChemists, nlambda=10) summary(fm_zip) fm_zinb <- zipath(art ~ . | 1, data = bioChemists, family = "negbin", nlambda=10) summary(fm_zinb) ## inflation with regressors ## ("art ~ . | ." is "art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment") fm_zip2 <- zipath(art ~ . | ., data = bioChemists, nlambda=10) summary(fm_zip2) fm_zinb2 <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) summary(fm_zinb2) ### non-penalized regression, compare with zeroinfl fm_zinb3 <- zipath(art ~ . | ., data = bioChemists, family = "negbin", lambda.count=0, lambda.zero=0, reltol=1e-12) summary(fm_zinb3) library("pscl") fm_zinb4 <- zeroinfl(art ~ . | ., data = bioChemists, dist = "negbin") summary(fm_zinb4) ### offset exposure <- rep(0.5, dim(bioChemists)[1]) fm_zinb <- zipath(art ~ . +offset(log(exposure))| ., data = bioChemists, family = "poisson", nlambda=10) coef <- coef(fm_zinb) ### offset can't be specified in predict function as it has been contained pred <- predict(fm_zinb) ## without inflation ## ("art ~ ." is "art ~ fem + mar + kid5 + phd + ment") fm_pois <- glmreg(art ~ ., data = bioChemists, family = "poisson") coef <- coef(fm_pois) fm_nb <- glmregNB(art ~ ., data = bioChemists) coef <- coef(fm_nb) ### high-dimensional #R CMD check --use-valgrind can be too time extensive for the following model #bioChemists <- cbind(matrix(rnorm(915*100), nrow=915), bioChemists) #fm_zinb <- zipath(art ~ . | ., data = bioChemists, family = "negbin", nlambda=10) ## End(Not run)
Fit zero-inflated regression models for count data via penalized maximum likelihood.
zipath_fit(X, Z, Y, weights, offsetx, offsetz, standardize=TRUE, intercept = TRUE, family = c("poisson", "negbin", "geometric"), link = c("logit", "probit", "cloglog", "cauchit", "log"), penalty = c("enet", "mnet", "snet"), start = NULL, y = TRUE, x = FALSE, nlambda=100, lambda.count=NULL, lambda.zero=NULL, type.path=c("active", "nonactive"), penalty.factor.count=NULL, penalty.factor.zero=NULL, lambda.count.min.ratio=.0001, lambda.zero.min.ratio=.1, alpha.count=1, alpha.zero=alpha.count, gamma.count=3, gamma.zero=gamma.count, rescale=FALSE, init.theta=NULL, theta.fixed=FALSE, EM=TRUE, maxit.em=200, convtype=c("count", "both"), maxit= 1000, maxit.theta =10, reltol = 1e-5, thresh=1e-6, eps.bino=1e-5, shortlist=FALSE, trace=FALSE, ...)zipath_fit(X, Z, Y, weights, offsetx, offsetz, standardize=TRUE, intercept = TRUE, family = c("poisson", "negbin", "geometric"), link = c("logit", "probit", "cloglog", "cauchit", "log"), penalty = c("enet", "mnet", "snet"), start = NULL, y = TRUE, x = FALSE, nlambda=100, lambda.count=NULL, lambda.zero=NULL, type.path=c("active", "nonactive"), penalty.factor.count=NULL, penalty.factor.zero=NULL, lambda.count.min.ratio=.0001, lambda.zero.min.ratio=.1, alpha.count=1, alpha.zero=alpha.count, gamma.count=3, gamma.zero=gamma.count, rescale=FALSE, init.theta=NULL, theta.fixed=FALSE, EM=TRUE, maxit.em=200, convtype=c("count", "both"), maxit= 1000, maxit.theta =10, reltol = 1e-5, thresh=1e-6, eps.bino=1e-5, shortlist=FALSE, trace=FALSE, ...)
X |
predictor matrix of the count model |
Z |
predictor matrix of the zero model |
Y |
response variable |
weights |
optional numeric vector of weights. |
offsetx |
optional numeric vector with an a priori known component to be included in the linear predictor of the count model. |
offsetz |
optional numeric vector with an a priori known component to be included in the linear predictor of the zero model. |
intercept |
Should intercept(s) be fitted (default=TRUE) or set to zero (FALSE) |
standardize |
Logical flag for x variable standardization, prior to
fitting the model sequence. The coefficients are always returned on
the original scale.
Default is |
family |
character specification of count model family (a log link is always used). |
link |
character specification of link function in the binary zero-inflation model (a binomial family is always used). |
y, x
|
logicals. If |
penalty |
penalty considered as one of |
start |
starting values for the parameters in the linear predictor. |
nlambda |
number of |
lambda.count |
A user supplied |
lambda.zero |
A user supplied |
type.path |
solution path with default value |
penalty.factor.count, penalty.factor.zero
|
These are numeric vectors with the same length as predictor variables. that multiply |
lambda.count.min.ratio, lambda.zero.min.ratio
|
Smallest value for |
alpha.count |
The elastic net mixing parameter for the count part of model. The default value 1 implies no L_2 penalty, as in LASSO. |
alpha.zero |
The elastic net mixing parameter for the zero part of model. The default value 1 implies no L_2 penalty, as in LASSO. |
gamma.count |
The tuning parameter of the |
gamma.zero |
The tuning parameter of the |
rescale |
logical value, if TRUE, adaptive rescaling |
init.theta |
The initial value of |
theta.fixed |
Logical value only used for |
EM |
Using |
convtype |
convergency type, default is for count component only for speedy computation |
maxit.em |
Maximum number of EM algorithm |
maxit |
Maximum number of coordinate descent algorithm |
maxit.theta |
Maximum number of iterations for estimating |
eps.bino |
a lower bound of probabilities to be claimed as zero, for computing weights and related values when |
reltol |
Convergence criteria, default value 1e-5 may be reduced to make more accurate yet slow |
thresh |
Convergence threshold for coordinate descent. Defaults value is |
shortlist |
logical value, if TRUE, limited results return |
trace |
If |
... |
Other arguments which can be passed to |
The algorithm fits penalized zero-inflated count data regression models using the coordinate descent algorithm within the EM algorithm.
The returned fitted model object is of class "zipath" and is similar
to fitted "glm" and "zeroinfl" objects. For elements such as "coefficients" a list is returned with elements for the zero and count component,
respectively.
If type.path="active", the algorithm iterates for a pair (lambda_count, lambda_zero) in a loop:
Step 1: For initial coefficients start_count of the count model and start_zero of the zero model, the EM algorithm is iterated until convergence for the active set with non-zero coefficients determined from start_count and start_zero, respectively.
Step 2: EM is iterated for all the predict variables once.
Step 3: If active set obtained from Step 2 is the same as in Step 1, stop; otherwise, repeat Step 1 and Step 2.
If type.path="nonactive", the EM algorithm iterates for a pair (lambda_count, lambda_zero) with all the predict variables until convergence.
A set of standard extractor functions for fitted model objects is available for
objects of class "zipath", including methods to the generic functions
print, coef,
logLik, residuals,
predict.
See predict.zipath for more details on all methods.
The program may terminate with the following message:
Error in: while (j <= maxit.em && !converged) { : Missing value, where TRUE/FALSE is necessaryCalls: zipathAdditionally: Warning:In glmreg_fit(Znew, probi, weights = weights, standardize = standardize, : saturated model, exiting ...Execution halted
One possible reason is that the fitted model is too complex for the data. There are two suggestions to overcome the error. One is to reduce the number of variables. Second, find out what lambda values caused the problem and omit them. Try with other lambda values instead.
An object of class "zipath", i.e., a list with components including
coefficients |
a list with elements |
residuals |
a vector of raw residuals (observed - fitted), |
fitted.values |
a vector of fitted means, |
weights |
the case weights used, |
terms |
a list with elements |
theta |
estimate of the additional |
loglik |
log-likelihood of the fitted model, |
family |
character string describing the count distribution used, |
link |
character string describing the link of the zero-inflation model, |
linkinv |
the inverse link function corresponding to |
converged |
logical value, TRUE indicating successful convergence of |
call |
the original function call |
formula |
the original formula |
levels |
levels of the categorical regressors |
model |
the full model frame (if |
y |
the response count vector (if |
x |
a list with elements |
Zhu Wang <[email protected]>
Zhu Wang, Shuangge Ma, Michael Zappitelli, Chirag Parikh, Ching-Yun Wang and Prasad Devarajan (2014) Penalized Count Data Regression with Application to Hospital Stay after Pediatric Cardiac Surgery, Statistical Methods in Medical Research. 2014 Apr 17. [Epub ahead of print]
Zhu Wang, Shuangge Ma, Ching-Yun Wang, Michael Zappitelli, Prasad Devarajan and Chirag R. Parikh (2014) EM for Regularized Zero Inflated Regression Models with Applications to Postoperative Morbidity after Cardiac Surgery in Children, Statistics in Medicine. 33(29):5192-208.
Zhu Wang, Shuangge Ma and Ching-Yun Wang (2015) Variable selection for zero-inflated and overdispersed data with application to health care demand in Germany, Biometrical Journal. 57(5):867-84.