RobustAFT-package(RobustAFT)
RobustAFT-package()所属R语言包：RobustAFT
Robust Accelerated Failure Time Model Fitting
强大的加速失效时间模型拟合
译者：生物统计家园网 机器人LoveR
描述----------Description----------
This package computes the truncated maximum likelihood regression estimates described in Marazzi and Yohai (2004) and 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. The main functions of this package are TML.noncensored and TML.censored.
这个软件包计算截断的最大似然回归估计在马拉斯和Yohai的（2004年）和Locatelli等人描述。 （2010年）。的误差分布假定约遵循高斯或log的Weibull分布。截止拒绝离群值是固定的或自适应。这个包的主要功能是TML.noncensored和TML.censored。
Details
详细信息----------Details----------
</table> This package depends on the libraries robustbase, survival, splines.
</ table>这个包所依赖的的库robustbase，生存，花键。
（作者）----------Author(s)----------
Alfio Marazzi <Alfio.Marazzi@chuv.ch>, Jean-Luc Muralti
Maintainer: A. Randriamiharisoa <Alex.Randriamiharisoa@chuv.ch>
参考文献----------References----------
Marazzi A., Yohai V. (2004). Adaptively truncated maximum likelihood regression with asymmetric errors. Journal of Statistical Planning and Inference, 122, 271-291.
Locatelli I., Marazzi A., Yohai V. (2010). Robust accelerated failure time regression. Computational Statistics and Data Analysis, 55, 874-887.
Marazzi A. (1993) Algorithm, Routines, and S functions for Robust Statistics. Wadsworth & Brooks/cole, Pacific Grove, California.
实例----------Examples----------
# Example 1. This is the example described in Marazzi and Yohai (2004).[实施例1。这是在马拉斯和Yohai的（2004）所描述的例子。]
# ---------[---------]
# The two following auxiliary functions, not included in the library, [下面的两个辅助功能，不包括在库中，]
#must be loaded.[必须被加载。]
# ------------- Auxiliary functions -------------------------------------[-------------辅助功能----------------------------------- - ]
SDmux.lw <- function(x,theta,sigma,COV){
# Standard deviation of the conditional mean estimate: log-Weibull case[条件均值估计的标准偏差：登录韦伯为例]
np <- length(theta); nc <- ncol(COV); nr <- nrow(COV)
if (np!=length(x)) cat("length(x) must be the same as length(theta)")
if (nc!=nr) cat("COV is not a square matrix")
if (nc!=(np+1)) cat("ncol(COV) must be the same as length(theta)+1")
log.mu.x <- t(x)%*%theta+lgamma(1+sigma) # log of conditional mean estimate[登录的条件均值估计]
mu.x <- exp(log.mu.x) # conditional mean estimate[条件均值估计]
dg <- digamma(1+sigma)
COV.TT <- COV[1:np,1:np]
Var.S <- COV[(np+1),(np+1)]
COV.TS <- COV[1:np,(np+1)]
V.mu.x <- t(x)%*%COV.TT%*%x + dg^2*Var.S + 2*dg*(t(x)%*%COV.TS)
SD.mu.x <- as.numeric((sqrt(V.mu.x))*mu.x)
SD.mu.x}
plt <- function(LOS,Cost,Adm,theta.fr,sigma.fr,sd0.fr,sd1.fr,theta.ml,
sigma.ml,sd0.ml,sd1.ml){
# Plot of the conditional mean and confidence intervals: log-Weibull case[图的条件均值和置信区间：log韦伯为例]
par(mfrow=c(1,2),oma=c(0,0,2,0))
plot(LOS,Cost,type="n")
points(LOS[Adm==0],Cost[Adm==0],pch=1)
points(LOS[Adm==1],Cost[Adm==1],pch=16,col=2)
x0t <- x0%*%theta.fr; x1t <- x1t <- x1%*%theta.fr
lines(l0,exp(x0t)*gamma(1+sigma.fr))
lines(l0,exp(x1t)*gamma(1+sigma.fr),col=2)
z0min <- exp(x0t)*gamma(1+sigma.fr)-2.576*sd0.fr
z0max <- exp(x0t)*gamma(1+sigma.fr)+2.576*sd0.fr
z1min <- exp(x1t)*gamma(1+sigma.fr)-2.576*sd1.fr
z1max <- exp(x1t)*gamma(1+sigma.fr)+2.576*sd1.fr
lines(l0,z0min,lty=2,col=1)
lines(l0,z0max,lty=2,col=1)
lines(l0,z1min,lty=2,col=1)
lines(l0,z1max,lty=2,col=1)
polygon(c(l0,rev(l0)), c(z0min,rev(z0max)), border=FALSE, density=10, angle=90)
polygon(c(l0,rev(l0)), c(z1min,rev(z1max)), border=FALSE, density=12, angle=90,col=2)
plot(LOS,Cost,type="n")
points(LOS[Adm==0],Cost[Adm==0],pch=1)
points(LOS[Adm==1],Cost[Adm==1],pch=16,col=2)
x0t <- x0%*%theta.ml; x1t <- x1t <- x1%*%theta.ml
lines(l0,exp(x0t)*gamma(1+sigma.ml))
lines(l0,exp(x1t)*gamma(1+sigma.ml),col=2)
z0min <- exp(x0t)*gamma(1+sigma.ml)-2.576*sd0.ml
z0max <- exp(x0t)*gamma(1+sigma.ml)+2.576*sd0.ml
z1min <- exp(x1t)*gamma(1+sigma.ml)-2.576*sd1.ml
z1max <- exp(x1t)*gamma(1+sigma.ml)+2.576*sd1.ml
lines(l0,z0min,lty=2,col=1)
lines(l0,z0max,lty=2,col=1)
lines(l0,z1min,lty=2,col=1)
lines(l0,z1max,lty=2,col=1)
polygon(c(l0,rev(l0)), c(z0min,rev(z0max)), border=FALSE, density=10, angle=90)
polygon(c(l0,rev(l0)), c(z1min,rev(z1max)), border=FALSE, density=12, angle=90,col=2)}
#------ End of auxiliary functions --------------------------------------------------[------结束的辅助功能---------------------------------------- ----------]
library(RobustAFT)
data(D243)
Cost <- D243$Cost # Cost (Swiss francs)[费用（瑞士法郎）]
LOS <- D243$LOS # Length of stay (days)[住院天数（天）]
Adm <- D243$Typadm; Adm <- (Adm==" Urg")*1 # Type of admission [入院类型]
# (0=on notification, 1=Emergency)[（0 =通知，紧急）]
Ass <- D243$Typass; Ass <- (Ass=="P" )*1 # Type of insurance (0=usual, 1=private)[保险类型（0 =往常一样，1 =私人）]
Age <- D243$age # Age (years)[年龄（岁）]
Dst <- D243$dest; Dst <- (Dst=="DOMI")*1 # Destination (1=Home, 0=another hospital)[目标（1 =首页，0 =另一家医院）]
Sex <- D243$Sexe; Sex <- (Sex=="M" )*1 # Sex (1=Male, 0=Female)[性别（1 = 0 =男，女）]
# Plot data[绘图数据]
par(mfrow=c(1,2))
plot(LOS,Cost); plot(log(LOS),log(Cost))
# log-Weibull fits[log - 韦伯适合]
# ----------------[----------------]
# Full robust model[完整的强大的模型]
zwff <- TML.noncensored(log(Cost)~log(LOS)+Adm+Ass+Age+Sex+Dst,
errors="logWeibull")
summary(zwff)
# Reduced model[减少模型]
zwfr <- update(zwff,log(Cost)~log(LOS)+Adm)
summary(zwfr)
# Residual plots[残差图]
par(mfrow=c(1,2))
plot(zwfr,which=c(1,3))
# Plot robust predictions on log-log scale[图强大的预测数坐标]
par(mfrow=c(1,1))
l0 <- seq(from=2,to=60,by=0.5)
x0 <- as.matrix(cbind(1,log(l0),0))
x1 <- as.matrix(cbind(1,log(l0),1))
plot(log(LOS),log(Cost),type="n")
points(log(LOS[Adm==1]),log(Cost[Adm==1]),pch=16,col=2)
points(log(LOS[Adm==0]),log(Cost[Adm==0]),pch=1)
lines(log(l0),predict(zwfr,x0))
lines(log(l0),predict(zwfr,x1),col=2)
# Maximum likelihood : full model[最大的可能性：完整的模型]
zmlf <- TML.noncensored(log(Cost)~log(LOS)+Adm+Ass+Age+Sex+Dst,
errors="logWeibull",cu=100)
summary(zmlf)
# Maximum likelihood : reduced model[最大的可能性：简化模型]
zmlr <- update(zmlf,log(Cost)~log(LOS)+Adm)
summary(zmlr)
# Plot conditional means and confidence intervals[图有条件的平均值和置信区间]
l0 <- seq(from=2,to=62,by=0.5)
x0 <- as.matrix(cbind(1,log(l0),0))
x1 <- as.matrix(cbind(1,log(l0),1))
theta.fr <- coef(zwfr)
sigma.fr <- zwfr$v1
COV.fr <- vcov(zwfr)
sd0.fr <- apply(x0,1,SDmux.lw,theta.fr,sigma.fr,COV.fr)
sd1.fr <- apply(x1,1,SDmux.lw,theta.fr,sigma.fr,COV.fr)
theta.ml <- coef(zmlr)
sigma.ml <- zmlr$v1
COV.ml <- zmlr$COV
sd0.ml <- apply(x0,1,SDmux.lw,theta.ml,sigma.ml,COV.ml)
sd1.ml <- apply(x1,1,SDmux.lw,theta.ml,sigma.ml,COV.ml)
plt(LOS,Cost,Adm,theta.fr,sigma.fr,sd0.fr,sd1.fr,theta.ml,sigma.ml,sd0.ml,sd1.ml)
# Gaussian fits (for comparison)[高斯适合（进行比较）]
# ------------------------------[------------------------------]
# Reduced model[减少模型]
zgfr <- TML.noncensored(log(Cost)~log(LOS)+Adm,errors="Gaussian")
summary(zgfr)
# Residual plots[残差图]
par(mfrow=c(1,2))
plot(zgfr,which=c(1,3))
# Classical Gaussian fit[古典高斯拟合]
lr <- lm(log(Cost)~log(LOS)+Adm)
summary(lr)
# Compare several fits[比较了几种适合]
# --------------------[--------------------]
comp <- fits.compare(TML.logWeibull=zwfr,ML.logWeibull=zmlr,least.squares=lr)
comp
plot(comp,leg.position=c("topleft","topleft","bottomleft")) # click on graphics[点击图形]
#[]
#------------------------------------------------------------------------------[-------------------------------------------------- ----------------------------]
#[]
# Example 2. This is the example described in Locatelli Marazzi and Yohai (2010).[实施例2。是洛卡特利马拉斯和Yohai的（2010）。]
# ---------[---------]
# This is the example described in Locatelli et al. (2010). [这是洛卡泰利等人所描述的例子。 （2010年）。]
# The estimates are slighty different due to changes in the algorithm for the [的估计略为不同的算法变化造成的]
# final estimate.[最终估计。]
# Remove data of Example 1[删除数据的实施例1中]
rm(Cost,LOS,Adm,Ass,Age,Dst,Sex)
data(MCI)
attach(MCI)
# Exploratory Analysis[探索性分析]
par(mfrow=c(1,1))
plot(Age,log(LOS),type= "n",cex=0.7)
# (1) filled square : regular, complete [（1）填方：定期，完整的]
# (2) empty square : regular, censored[（2）空方：定期，删]
# (3) filled triangle : emergency, complete[（3）填充三角形的紧急情况下，]
# (4) empty triangle : emergency, censored[（4）空三角：紧急情况下，审查]
points(Age[Dest==1 & TypAdm==0], log(LOS)[Dest==1 & TypAdm==0], pch=15,cex=0.7) # (1)[（1）]
points(Age[Dest==0 & TypAdm==0], log(LOS)[Dest==0 & TypAdm==0], pch=0, cex=0.7) # (2)[（2）]
points(Age[Dest==1 & TypAdm==1], log(LOS)[Dest==1 & TypAdm==1], pch=17,cex=0.7) # (3)[（3）]
points(Age[Dest==0 & TypAdm==1], log(LOS)[Dest==0 & TypAdm==1], pch=2, cex=0.7) # (4)[（4）]
# Maximum Likelihood[最大似然法（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")
#[]
detach(MCI)
转载请注明:出自 生物统计家园网(http://www.biostatistic.net)。
注：
注1：为了方便大家学习，本文档为生物统计家园网机器人LoveR翻译而成，仅供个人R语言学习参考使用，生物统计家园保留版权。
注2：由于是机器人自动翻译，难免有不准确之处，使用时仔细对照中、英文内容进行反复理解，可以帮助R语言的学习。
注3：如遇到不准确之处，请在本贴的后面进行回帖，我们会逐渐进行修订。
|
发表于 2012-9-27 22:04:19
