2017-07-13 48 views
1

我不能強制rc樣條函數預測參考在函數內部改變。我希望65歲是參考(yhat = 1,lower = 1,upper = 1)。當代碼在「分析」功能之外時,它可以很好地工作。即使dd $限制確實更改爲參考65(請參閱輸出),我懷疑它在函數內部時不會發生「更新」。可能不同的環境?我嘗試了幾個小時來設置環境等,不幸的是沒有成功。任何幫助將不勝感激!在r中改變參考值限制三次樣條函數

library(Hmisc) 
    library(survival) 
    library(rms) 
    library(cmprsk) 

    bk.tst <- function(analysis) 
    { 
    if (analysis=="test") { 
     dt<-m 
    } 
    outcomes<-c("compos", "acs", "death", "vascdeath", "stroke", "majbleed") 
    for (i in outcomes) { 
    fvl<- dt$fvl 
    age<-dt$age 
    ndt <- data.frame(age, fvl) 
    dd<-assign('dd', datadist(ndt), pos=1)  
    options(datadist='dd') 
    tm<-paste("tt",i, sep="") 
    SurvObj <- with(dt, Surv(eval(parse(text=tm)), eval(parse(text=i))==1)) 
    f<-cph(SurvObj ~ fvl*rcs(age,c(60,70,80)), type="Survival", method="exact", x=T, y=T) 
    print(dd$limits) 
    dd$limits["Adjust to","age"] <- 65 
    print(dd$limits) 
    g <- update(f) 
    ano<-anova(f) 
    age.intr.rcsplines<-Predict(g, age=45:85, fvl, ref.zero=TRUE, fun=exp) 
    print(age.intr.rcsplines[20:22,]) 
    } 
    } 

    bk.tst("test") 

#the output is: 
        age fvl 
Low:effect  57.01027 0 
Adjust to  64.62697 0 
High:effect  71.69884 1 
Low:prediction 39.48545 0 
High:prediction 84.64122 1 
Low    35.80287 0 
High   92.15606 1 
        age fvl 
Low:effect  57.01027 0 
Adjust to  65.00000 0 
High:effect  71.69884 1 
Low:prediction 39.48545 0 
High:prediction 84.64122 1 
Low    35.80287 0 
High   92.15606 1 
    age fvl  yhat  lower upper 
20 64 0 0.986740 0.9599284 1.014300 
21 65 0 1.008632 0.9929697 1.024541 
22 66 0 1.034987 0.9799775 1.093084 

# when the code is outside the function analysis, then the output is correct with 65 as a reference. That is what I want to happen also inside the function "analysis": 
        age fvl 
Low:effect  57.01027 0 
Adjust to  64.62697 0 
High:effect  71.69884 1 
Low:prediction 39.48545 0 
High:prediction 84.64122 1 
Low    35.80287 0 
High   92.15606 1 
        age fvl 
Low:effect  57.01027 0 
Adjust to  65.00000 0 
High:effect  71.69884 1 
Low:prediction 39.48545 0 
High:prediction 84.64122 1 
Low    35.80287 0 
High   92.15606 1 
    age fvl  yhat  lower upper 
20 64 0 0.9782956 0.9369436 1.021473 
21 65 0 1.0000000 1.0000000 1.000000 
22 66 0 1.0261294 0.9868852 1.066934 

回答

0

最後,在函數bk.tst中定義環境如下所示工作!

bk.tst <- function(analysis, env = parent.frame()) 
{ 
if (analysis=="test") { 
    dt<-m 
} 
outcomes<-c("compos", "acs", "death", "vascdeath", "stroke", "majbleed") 
for (i in outcomes) { 
fvl<- dt$fvl 
age<-dt$age 
env$ndt <- data.frame(age, fvl) 
env$dd<-datadist(env$ndt)  
options(datadist='dd') 
tm<-paste("tt",i, sep="") 
SurvObj <- with(dt, Surv(eval(parse(text=tm)), eval(parse(text=i))==1)) 
env$f<-cph(SurvObj ~ fvl*rcs(age,c(60,70,80)), type="Survival", method="exact", x=T, y=T) 
print(env$dd$limits) 
env$dd$limits["Adjust to","age"] <- 65 
print(env$dd$limits) 
g <- update(env$f) 
ano<-anova(g) 
age.intr.rcsplines<-Predict(g, age=45:85, fvl, ref.zero=TRUE, fun=exp) 
print(age.intr.rcsplines[20:22,]) 
} 
} 

bk.tst("test")