R: Truncated Maximum Likelihood Regression With Censored...
TML.censored
R Documentation
Truncated Maximum Likelihood Regression With Censored Observations
Description
This function computes the truncated maximum likelihood estimates of accelerated failure time regression
described in Locatelli et al. (2010).
The error distribution is assumed to follow approximately a Gaussian or a log-Weibull distribution.
The cut-off values for outlier rejection are fixed or adaptive.
A formula, i.e., a symbolic description of the model
to be adjusted (cf. glm or lm).
data
An optional data frame containing the variables in the model. If not
found in data, the variables are taken from environment(formula),
typically the environment from which robaft is called.
delta
Vector of 0 and 1.
0: censored observation.
1: complete observation.
errors
"Gaussian": the error distribution is assumed to be Gaussian.
"logWeibull" : the error distribution is assumed to be log-Weibull.
initial
"S": initial S-estimate.
"input": the initial estimate is given on input.
input
A list(theta=c(...),sigma=...): initial input estimates where
theta is a vector of p coefficients and sigma a scalar scale.
Required when initial="input".
otp
"adaptive": adaptive cut-off.
"fixed": non adaptive cut-off.
cov
If TRUE the covariance matrix is computed.
cu
Preliminary minimal upper cut-off.
The default is 2.5 in the Gaussian case and 1.855356 in the log-Weibull case.
control.S
A list of control parameters for the computation of the initial S estimates.
See the function TML.censored.control.S for the default values.
control.ref
A list of control parameters for the refinement algorithm of the initial S estimates.
See the function TML.censored.control.ref for the default values.
control.tml
AA list of control parameters for the computation of the final estimates.
See the function TML.censored.control.tml for the default values.
Value
TML.censored returns an object of class "TML".
The function summary can be used to obtain or print a summary of the results.
The generic extractor functions fitted, residuals and
weights can be used to extract various elements of the object returned
by TML.censored. The function update can be used to update the model.
An object of class "TML" is a list with at least the following components:
th0
Initial coefficient estimates.
v0
Initial scale estimate.
nit.ref
Reached number of iteration in the refinement step for the initial estimates.
th1
Final coefficient estimates.
v1
Final scale estimate.
nit.tml
Number of iterations reached in IRLS algorithm for the final estimates.
tu,tl
Final cut-off values.
alpha
Estimated proportion of retained observations.
tn
Number of retained observations.
weights
Vector of weights (0 for rejected observations, 1 for retained observations).
COV
Covariance matrix of the final estimates (th1[1],...,th1[p],v1) (where p=ncol(X)).
residuals
Residuals of noncensored observations are calculated as response minus fitted values.
For censored observations, the the expected residuals given that
the response is larger than the recorded censored value are provided.
fitted.values
The fitted mean values.
call
The matched call.
formula
The formula supplied.
terms
The terms object used.
data
The data argument.
References
Locatelli I., Marazzi A., Yohai V. (2010). Robust accelerated failure time regression.
Computational Statistics and Data Analysis, 55, 874-887.
# This is the example described in Locatelli et al. (2010).
# The estimates are slighty different than those of the paper due to changes
# in the algorithm for the final estimate.
#
## Not run:
data(MCI)
attach(MCI)
# Exploratory Analysis
plot(Age,log(LOS),type= "n",cex=0.7)
# (1) filled square : regular, complete
# (2) empty square : regular, censored
# (3) filled triangle : emergency, complete
# (4) empty triangle : emergency, censored
points(Age[Dest==1 & TypAdm==0], log(LOS)[Dest==1 & TypAdm==0], pch=15,cex=0.7) # (1)
points(Age[Dest==0 & TypAdm==0], log(LOS)[Dest==0 & TypAdm==0], pch=0, cex=0.7) # (2)
points(Age[Dest==1 & TypAdm==1], log(LOS)[Dest==1 & TypAdm==1], pch=17,cex=0.7) # (3)
points(Age[Dest==0 & TypAdm==1], log(LOS)[Dest==0 & TypAdm==1], pch=2, cex=0.7) # (4)
# Maximum Likelihood
ML <- survreg(Surv(log(LOS), Dest) ~ TypAdm*Age, dist="gaussian")
summary(ML)
B.ML <- ML$coef
S.ML <- ML$scale
abline(c(B.ML[1] ,B.ML[3] ),lwd=1,col="grey",lty=1)
abline(c(B.ML[1]+B.ML[2],B.ML[3]+B.ML[4]),lwd=1,col="grey",lty=1)
# Robust Accelerated Failure Time Regression with Gaussian errors
ctrol.S <- list(N=150, q=5, sigma0=1, MAXIT=100, TOL=0.001,seed=123)
ctrol.ref <- list(maxit.sigma=2,tol.sigma=0.0001,maxit.Beta=2,tol.Beta=0.0001,
Maxit.S=50, tol.S.sigma=0.001, tol.S.Beta=0.001,alg.sigma=1,nitmon=FALSE)
ctrol.tml <- list(maxit.sigma=50,tol.sigma=0.0001,maxit.Beta=50,tol.Beta=0.0001,
Maxit.TML=50, tol.TML.sigma=0.001, tol.TML.Beta=0.001, alg.sigma=1,nitmon=FALSE)
WML<-TML.censored(log(LOS)~TypAdm*Age,data=MCI,delta=Dest,otp="adaptive",
control.S=ctrol.S,control.ref=ctrol.ref,control.tml=ctrol.tml)
summary(WML)
B.WML<-coef(WML)
abline(c(B.WML[1] ,B.WML[3] ),lty=1, col="red")
abline(c(B.WML[1]+B.WML[2],B.WML[3]+B.WML[4]),lty=1, col="red")
## End(Not run)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(RobustAFT)
Loading required package: robustbase
Loading required package: survival
Attaching package: 'survival'
The following object is masked from 'package:robustbase':
heart
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/RobustAFT/TML.censored.rd_%03d_medium.png", width=480, height=480)
> ### Name: TML.censored
> ### Title: Truncated Maximum Likelihood Regression With Censored
> ### Observations
> ### Aliases: TML.censored
> ### Keywords: Regression Robust Accelerated Failure Time
>
> ### ** Examples
>
> # This is the example described in Locatelli et al. (2010).
> # The estimates are slighty different than those of the paper due to changes
> # in the algorithm for the final estimate.
> #
> ## Not run:
> ##D data(MCI)
> ##D attach(MCI)
> ##D
> ##D # Exploratory Analysis
> ##D plot(Age,log(LOS),type= "n",cex=0.7)
> ##D
> ##D # (1) filled square : regular, complete
> ##D # (2) empty square : regular, censored
> ##D # (3) filled triangle : emergency, complete
> ##D # (4) empty triangle : emergency, censored
> ##D
> ##D points(Age[Dest==1 & TypAdm==0], log(LOS)[Dest==1 & TypAdm==0], pch=15,cex=0.7) # (1)
> ##D points(Age[Dest==0 & TypAdm==0], log(LOS)[Dest==0 & TypAdm==0], pch=0, cex=0.7) # (2)
> ##D points(Age[Dest==1 & TypAdm==1], log(LOS)[Dest==1 & TypAdm==1], pch=17,cex=0.7) # (3)
> ##D points(Age[Dest==0 & TypAdm==1], log(LOS)[Dest==0 & TypAdm==1], pch=2, cex=0.7) # (4)
> ##D
> ##D # Maximum Likelihood
> ##D ML <- survreg(Surv(log(LOS), Dest) ~ TypAdm*Age, dist="gaussian")
> ##D summary(ML)
> ##D B.ML <- ML$coef
> ##D S.ML <- ML$scale
> ##D
> ##D abline(c(B.ML[1] ,B.ML[3] ),lwd=1,col="grey",lty=1)
> ##D abline(c(B.ML[1]+B.ML[2],B.ML[3]+B.ML[4]),lwd=1,col="grey",lty=1)
> ##D
> ##D # Robust Accelerated Failure Time Regression with Gaussian errors
> ##D ctrol.S <- list(N=150, q=5, sigma0=1, MAXIT=100, TOL=0.001,seed=123)
> ##D
> ##D ctrol.ref <- list(maxit.sigma=2,tol.sigma=0.0001,maxit.Beta=2,tol.Beta=0.0001,
> ##D Maxit.S=50, tol.S.sigma=0.001, tol.S.Beta=0.001,alg.sigma=1,nitmon=FALSE)
> ##D
> ##D ctrol.tml <- list(maxit.sigma=50,tol.sigma=0.0001,maxit.Beta=50,tol.Beta=0.0001,
> ##D Maxit.TML=50, tol.TML.sigma=0.001, tol.TML.Beta=0.001, alg.sigma=1,nitmon=FALSE)
> ##D
> ##D WML<-TML.censored(log(LOS)~TypAdm*Age,data=MCI,delta=Dest,otp="adaptive",
> ##D control.S=ctrol.S,control.ref=ctrol.ref,control.tml=ctrol.tml)
> ##D
> ##D summary(WML)
> ##D
> ##D B.WML<-coef(WML)
> ##D abline(c(B.WML[1] ,B.WML[3] ),lty=1, col="red")
> ##D abline(c(B.WML[1]+B.WML[2],B.WML[3]+B.WML[4]),lty=1, col="red")
> ## End(Not run)
>
>
>
>
>
> dev.off()
null device
1
>