BUPA Liver Disorder Data

Description

345 male patients' blood test result and liver disorder status.

Usage

data(BUPA)

Format

Details

All the variables, except for the response, have been scaled to [0,1] interval. To transform back to the original scale, use the formula:

x=min+(max-min)*z.

Source

Richard S. Forsyth at BUPA Medical Research Ltd.

Compute adaptive weights by fitting a SS-ANOVA model

Description

A preliminary estimate \tilde{\eta} is first obtained by fitting a smoothing spline ANOVA model, and then use the inverse L_2-norm, ||\tilde{\eta}_j||^{-\gamma}, as the initial weight for the j-th functional component.

Usage

 SSANOVAwt(x,y,tau,family=c("Gaussian","Binomial","Cox","Quantile"),mscale=rep(1,ncol(x)),
               gamma=1,scale=FALSE,nbasis,basis.id,cpus) 

Arguments

Details

The initial mean function is estimated via a smooothing spline objective function. In the SS-ANOVA model framework, the regression function is assumed to have an additive form

\eta(x)=b+\sum_{j=1}^p\eta_j(x^{(j)}),

where b denotes intercept and \eta_j denotes the main effect of the j-th covariate.

For "Gaussian" response, the mean regression function is estimated by minimizing the objective function:

\sum_i(y_i-\eta(x_i))^2/nobs+\lambda_0\sum_{j=1}^p\alpha_j||\eta_j||^2.

where RSS is residual sum of squares.

For "Binomial" response, the regression function is estimated by minimizing the objective function:

-log-likelihood/nobs+\lambda_0\sum_{j=1}^p\alpha_j||\eta_j||^2

For "Quantile" regression model, the quantile function, is estimated by minimizing the objective function:

\sum_i\rho(y_i-\eta(x_i))/nobs+\lambda_0\sum_{j=1}^p\alpha_j||\eta_j||^2.

For "Cox" regression model, the log-hazard function, is estimated by minimizing the objective function:

-log-Partial Likelihood/nobs+\lambda_0\sum_{j=1}^p\alpha_j||\eta_j||^2.

The smoothing parameter \lambda_0 is tuned by 5-fold Cross-Validation, if family="Gaussian", "Binomial" or "Quantile", and Approximate Cross-Validation, if family="Cox". But the smoothing parameters \alpha_j are given in the argument mscale.

The adaptive weights are then fiven by ||\tilde{\eta}_j||^{-\gamma}.

Value

Author(s)

Hao Helen Zhang and Chen-Yen Lin

References

Storlie, C. B., Bondell, H. D., Reich, B. J. and Zhang, H. H. (2011) "Surface Estimation, Variable Selection, and the Nonparametric Oracle Property", Statistica Sinica, 21, 679–705.

Examples

## Adaptive COSSO Model
## Binomial
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*7,0,1),nc=7))
trueProb=1/(1+exp(-x[,1]-sin(2*pi*x[,2])-5*(x[,4]-0.4)^2))
y=rbinom(200,1,trueProb)

Binomial.wt=SSANOVAwt(x,y,family="Bin")
ada.B.Obj=cosso(x,y,wt=Binomial.wt,family="Bin")

## Not run: 
## Gaussian
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*7,0,1),nc=7))
y=x[,1]+sin(2*pi*x[,2])+5*(x[,4]-0.4)^2+rnorm(200,0,1)
Gaussian.wt=SSANOVAwt(designx,response,family="Gau")
ada.G.Obj=cosso(x,y,wt=Gaussian.wt,family="Gaussian")

## End(Not run)

Fit a generalized nonparametric model with cosso penalty

Description

A comprehensive method for fitting various type of regularized nonparametric regression models using cosso penalty. Fits mean, logistic, Cox and quantile regression.

Usage

cosso(x,y,tau,family=c("Gaussian","Binomial","Cox","Quantile"),wt=rep(1,ncol(x)),
      scale=FALSE,nbasis,basis.id,cpus)

Arguments

Details

In the SS-ANOVA model framework, the regression function is assumed to have an additive form

\eta(x)=b+\sum_{j=1}^p\eta_j(x^{(j)}),

where b denotes intercept and \eta_j denotes the main effect of the j-th covariate.

For "Gaussian" response, the mean function is estimated by minimizing the objective function:

\sum_i(y_i-\eta(x_i))^2/nobs+\lambda_0\sum_{j=1}^p\theta^{-1}_jw_j^2||\eta_j||^2, s.t.~\sum_{j=1}^p\theta_j\leq M.

For "Binomial" response, the log-odd function is estimated by minimizing the objective function:

-log-likelihood/nobs+\lambda_0\sum_{j=1}^p\theta^{-1}_jw_j^2||\eta_j||^2, s.t.~\sum_{j=1}^p\theta_j\leq M.

For "Quantile" regression model, the quantile function, is estimated by minimizing the objective function:

\sum_i\rho_{\tau}(y_i-\eta(x_i))/nobs+\lambda_0\sum_{j=1}^p\theta^{-1}_jw_j^2||\eta_j||^2, s.t.~\sum_{j=1}^p\theta_j\leq M.

For "Cox" regression model, the log-relative hazard function is estimated by minimizing the objective function:

-log-Partial Likelihood/nobs+\lambda_0\sum_{j=1}^p\theta^{-1}_jw_j^2||\eta_j||^2, s.t.~\sum_{j=1}^p\theta_j\leq M.

For identifiability sake, the intercept term in Cox model is absorbed into basline hazard, or equivalently set b=0.

For large data sets, we can reduce the computational load of the optimization problem by selecting a subset of the observations of size nbais as "knots", which reduces the dimension of the kernel matrices from nobs to nbasis. Unless specified via basis.id or nbasis, the default number of "knots" is max(40,12*nobs^(2/9)) for "Gaussian" and "Binomial" and max(35,11 * nobs^(2/9)) for "Cox".

The weights can be specified based on either user's own discretion or adaptively computed from initial function estimates. See Storlie et al. (2011) for more discussions. One possible choice is to specify the weights as the inverse L_2 norm of initial function estimator, see SSANOVAwt.

Value

An object with S3 class "cosso".

Author(s)

Hao Helen Zhang and Chen-Yen Lin

References

Lin, Y. and Zhang, H. H. (2006) "Component Selection and Smoothing in Smoothing Spline Analysis of Variance Models", Annals of Statistics, 34, 2272–2297.

Leng, C. and Zhang, H. H. (2006) "Model selection in nonparametric hazard regression", Nonparametric Statistics, 18, 417–429.

Zhang, H. H. and Lin, Y. (2006) "Component Selection and Smoothing for Nonparametric Regression in Exponential Families", Statistica Sinica, 16, 1021–1041.

Storlie, C. B., Bondell, H. D., Reich, B. J. and Zhang, H. H. (2011) "Surface Estimation, Variable Selection, and the Nonparametric Oracle Property", Statistica Sinica, 21, 679–705.

See Also

plot.cosso, predict.cosso, tune.cosso

Examples

## Gaussian
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*9,0,1),nc=9))
y=x[,1]+sin(2*pi*x[,2])+5*(x[,4]-0.4)^2+rnorm(200,0,1)
G.Obj=cosso(x,y,family="Gaussian")
plot.cosso(G.Obj,plottype="Path")

## Not run: 
## Use all observations as knots
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*9,0,1),nc=9))
y=x[,1]+sin(2*pi*x[,2])+5*(x[,4]-0.4)^2+rnorm(200,0,1)
G.Obj=cosso(x,y,family="Gaussian",nbasis=200)
## Clean up
rm(list=ls())

## Binomial
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*9,0,1),nc=9))
trueProb=1/(1+exp(-x[,1]-sin(2*pi*x[,2])-5*(x[,4]-0.4)^2))
y=rbinom(200,1,trueProb)

B.Obj=cosso(x,y,family="Bin")
## Clean up
rm(list=ls())

## Cox
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*9,0,1),nc=9))
hazard=x[,1]+sin(2*pi*x[,2])+5*(x[,4]-0.4)^2
surTime=rexp(200,exp(hazard))
cenTime=rexp(200,exp(-hazard)*runif(1,4,6))
y=cbind(time=apply(cbind(surTime,cenTime),1,min),status=1*(surTime<cenTime))
C.obj=cosso(x,y,family="Cox",cpus=1)

## Try parallel computing
C.obj=cosso(x,y,family="Cox",cpus=4)

## Clean up
rm(list=ls())

## Quantile
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*7,0,1),nc=7))
y=x[,1]+sin(2*pi*x[,2])+5*(x[,4]-0.4)^2+rt(200,3)
Q.obj=cosso(x,y,0.3,family="Quan",cpus=1)

## Try parallel computing
Q.obj=cosso(x,y,0.3,family="Quan",cpus=4)

## End(Not run)

Internal Functions

Description

Internal cosso functions.

Details

These are not to be called by the user (or in some cases are just waiting for proper documentation to be written :).

Ozone pollution data in Los Angels, 1976

Description

This is the ozone data used in Breiman and Friedman (1985). This dataset contains 330 observations, and each observation is a daily measurement.

Usage

data(ozone)

Format

Details

All the variables, except for the response, have been scaled to [0,1] interval. To transform back to the original scale, use the formula:

x=min+(max-min)*z.

Source

Breiman, L. and Friedman, J. (1985), "Estimating Optimal Transformations for Multiple Regression and Correlation", Journal of the American Statistical Association, 80, 580–598.

Plot method for COSSO object

Description

Plot L_2 norm solution path or main effects of selected functional components

Usage

## S3 method for class 'cosso'
plot(x,M,plottype =c("Path","Functionals"),eps=1e-7,...)

Arguments

Value

NULL

Author(s)

Hao Helen Zhang and Chen-Yen Lin

See Also

predict.cosso

Examples

set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*7,0,1),nc=7))
trueProb=1/(1+exp(-x[,1]-sin(2*pi*x[,2])-5*(x[,4]-0.4)^2))
y=rbinom(200,1,trueProb)

B.Obj=cosso(x,y,family="Bin")
plot.cosso(B.Obj,plottype="Path")
plot.cosso(B.Obj,M=2,plottype="Func")

Make predictions or extract coefficients from a cosso model

Description

Make prediction for future observations or extract the model parameters at a particular smoothing parameter.

Usage

 
## S3 method for class 'cosso'
predict(object,xnew,M,type=c("fit","coefficients","nonzero"),eps=1e-7,...) 

Arguments

Value

The object returned depends on type.

When type="fit", predicted

eta

function value will be given at the new design points xnew.

When type="coefficients", three sets of coefficients will be returned.

When type="nonzero", a list of the indices of the nonconstant functional components will be returned.

Author(s)

Hao Helen Zhang and Chen-Yen Lin

See Also

plot.cosso

Examples

 
## Gaussian
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*7,0,1),nc=7))
y=x[,1]+sin(2*pi*x[,2])+5*(x[,4]-0.4)^2+rnorm(200,0,1)
G.Obj=cosso(x,y,family="Gaussian")
predict.cosso(G.Obj,M=2,type="nonzero")
predict.cosso(G.Obj,xnew=x[1:3,],M=2,type="fit")
## Clean up
rm(list=ls())

## Not run: 
## Binomial
set.seed(20130310)
x=cbind(rbinom(200,1,.7),matrix(runif(200*9,0,1),nc=9))
trueProb=1/(1+exp(-x[,1]-sin(2*pi*x[,2])-5*(x[,4]-0.4)^2))
y=rbinom(200,1,trueProb)

B.Obj=cosso(x,y,family="Bin") 
f.hat=predict.cosso(B.Obj,xnew=x,M=2,type="fit")
prob.hat=1/(1+exp(-f.hat))
## Clean up
rm(list=ls())

## End(Not run)

Tuning procedure for cosso

Description

Compute K-fold cross-validated score and plot cross-validated score against a grid values of smooth parameter M.

Usage

 
tune.cosso(object,folds=5,plot.it=TRUE)   

Arguments

Value

Author(s)

Hao Helen Zhang and Chen-Yen Lin

See Also

cosso, predict.cosso

Examples

 
## Binomial
set.seed(20130310)
x=cbind(rbinom(150,1,.7),matrix(runif(150*5,0,1),nc=5))
trueProb=1/(1+exp(-x[,1]-sin(2*pi*x[,2])-5*(x[,4]-0.4)^2))
y=rbinom(150,1,trueProb)

B.Obj=cosso(x,y,family="Bin",nbasis=30) 
tune.cosso(B.Obj,4,TRUE)

Veterans' Administration Lung Cancer study

Description

Randomized trial of two treatment regimens for lung cancer.

Usage

data(veteran)

Format

Details

All the variables, except for the response, have been scaled to [0,1] interval. To transform back to the original scale, use the formula:

x=min+(max-min)*z.

Source

Kalbfleisch, J. and Prentice, R.L. (2002), The Statistical Analysis of Failure Time Data (Second Edition) Wiley: New Jersey.