2012-01-25 4 views
4

survfit()basehaz()을 함수 내에 사용하고 싶지만 작동하지 않습니다. 이 문제를 살펴볼 수 있습니까? 당신의 도움을 주셔서 감사합니다.함수 내부에서 수식 오류가 발생했습니다.

library(survival) 

n <- 50  # total sample size 
nclust <- 5 # number of clusters 
clusters <- rep(1:nclust,each=n/nclust) 
beta0 <- c(1,2) 
set.seed(13) 
#generate phmm data set 
Z <- cbind(Z1=sample(0:1,n,replace=TRUE), 
     Z2=sample(0:1,n,replace=TRUE), 
     Z3=sample(0:1,n,replace=TRUE)) 
b <- cbind(rep(rnorm(nclust),each=n/nclust),rep(rnorm(nclust),each=n/nclust)) 
Wb <- matrix(0,n,2) 
for(j in 1:2) Wb[,j] <- Z[,j]*b[,j] 
Wb <- apply(Wb,1,sum) 
T <- -log(runif(n,0,1))*exp(-Z[,c('Z1','Z2')]%*%beta0-Wb) 
C <- runif(n,0,1) 
time <- ifelse(T<C,T,C) 
event <- ifelse(T<=C,1,0) 
mean(event) 
phmmd <- data.frame(Z) 
phmmd$cluster <- clusters 
phmmd$time <- time 
phmmd$event <- event 

fmla <- as.formula("Surv(time, event) ~ Z1 + Z2") 

BaseFun <- function(x){ 
start.coxph <- coxph(x, phmmd) 

print(start.coxph) 

betahat <- start.coxph$coefficient 
print(betahat) 
print(333) 
print(survfit(start.coxph))                                                          
m <- basehaz(start.coxph) 
print(m) 
} 
BaseFun(fmla) 
Error in formula.default(object, env = baseenv()) : invalid formula 

을하지만, 다음과 같은 기능이 작동 : 다음 코드는 오류에 이르게

fit <- coxph(fmla, phmmd)  
basehaz(fit) 
+0

http://stackoverflow.com/q/10176524/210673 또한 동일한 문제로 보입니다. – Aaron

답변

5

그것은 범위 지정의 문제입니다. basehaz의 환경은 주의 사항 : 한편

environment(basehaz) 
<environment: namespace:survival> 

:

environment(BaseFun) 
<environment: R_GlobalEnv> 

따라서 그 기능 basehaz 함수 내부의 지역 변수를 찾을 수없는 이유입니다.

가능한 솔루션은 assign를 사용하여 상단에있는 X를 보내는 것입니다 :

BaseFun <- function(x){ 

    assign('x',x,pos=.GlobalEnv) 

    start.coxph <- coxph(x, phmmd) 
    print(start.coxph) 

    betahat <- start.coxph$coefficient 
    print(betahat) 
    print(333) 
    print(survfit(start.coxph)) 

    m <- basehaz(start.coxph) 
    print(m) 
    rm(x) 

     } 
    BaseFun(fmla) 

다른 솔루션보다 직접적으로 환경을 다루는 포함 할 수있다.

+0

좋은 찾으십시오. 지구 환경에 'm'을 부여 할 특별한 이유가 있습니까? –

+0

@ RomanLuštrik 실제로 오자 인 이유는 없습니다. 잡기위한 txs. – aatrujillob

+1

@AndresT 대단히 감사합니다. 그것은 작동합니다. Terry Therneau는 또한 해결책을 제시했습니다. 생존 루틴은 모델 행렬을 재구성해야하며 기본적으로 R에서 모델 수식이 처음 정의 된 컨텍스트에서 수행됩니다. 불행히도이 함수 밖에서 문제가 발생합니다. 인수 "x"는 외부 환경에서 알 수 없습니다. "해결 방법은 coxph 호출에"model = TRUE "를 추가하여 모델 프레임을 저장하고 survfit을 재구성 할 필요가 없도록하는 것입니다." – moli

2

나는 @ a1rujillob의 대답에 @ moli의 의견을 후속 조치 중입니다. 그것들은 도움이 되었기 때문에 그것이 나를 위해 어떻게 해결되었고, rpartpartykit 패키지와 비슷한 문제를 해결할 수있을 것이라고 생각했습니다.

일부 장난감 데이터 :

N <- 200 
data <- data.frame(X = rnorm(N),W = rbinom(N,1,0.5)) 
data <- within(data, expr = { 
    trtprob <- 0.4 + 0.08*X + 0.2*W -0.05*X*W 
    Trt <- rbinom(N, 1, trtprob) 
    outprob <- 0.55 + 0.03*X -0.1*W - 0.3*Trt 
    Outcome <- rbinom(N,1,outprob) 
    rm(outprob, trtprob) 
}) 

내가 훈련 (train_data) 및 테스트 세트에 데이터를 분할하고, train_data에 분류 트리를 양성하고자합니다.

다음은 사용하려는 수식과 다음 예제의 문제입니다. 이 수식을 정의하면 train_data 개체가 아직 존재하지 않습니다.

다음은 원래 기능과 유사한 내 기능입니다. 다시

badFunc <- function(data, my_formula){ 
    train_data <- data[1:100,] 
    ct_train <- rpart::rpart(
    data= train_data, 
    formula = my_formula, 
    method = "class") 
    ct_party <- partykit::as.party(ct_train) 
} 

이 함수를 실행하려고하면 OP와 비슷한 오류가 발생합니다.

library(rpart) 
library(partykit) 

bad_out <- badFunc(data=data, my_formula = my_formula) 
# Error in is.data.frame(data) : object 'train_data' not found 
# 10. is.data.frame(data) 
# 9. model.frame.default(formula = Trt ~ W + X, data = train_data, 
#   na.action = function (x) {Terms <- attr(x, "terms") ... 
# 8. stats::model.frame(formula = Trt ~ W + X, data = train_data, 
#   na.action = function (x) {Terms <- attr(x, "terms") ... 
# 7. eval(expr, envir, enclos) 
# 6. eval(mf, env) 
# 5. model.frame.rpart(obj) 
# 4. model.frame(obj) 
# 3. as.party.rpart(ct_train) 
# 2. partykit::as.party(ct_train) 
# 1. badFunc(data = data, my_formula = my_formula) 

print(bad_out) 
# Error in print(bad_out) : object 'bad_out' not found 

다행히 rpart()는 이러한 문제를 해결하기 위해 인수 model=TRUE을 지정할 수 있다는 점에서 coxph() 같다. 여기에 또 다른 논쟁이 있습니다.

goodFunc <- function(data, my_formula){ 
    train_data <- data[1:100,] 
    ct_train <- rpart::rpart(
    data= train_data, 
    ## This solved it for me 
    model=TRUE, 
    ## 
    formula = my_formula, 
    method = "class") 
    ct_party <- partykit::as.party(ct_train) 
} 
good_out <- goodFunc(data=data, my_formula = my_formula) 
print(good_out)  
# Model formula: 
# Trt ~ W + X 
# 
# Fitted party: 
# [1] root 
# | [2] X >= 1.59791: 0.143 (n = 7, err = 0.9) 
##### etc 

문서 rpart()에서 model 인수 : 그들은 (나에게) 항상 자연이 아닌 방법으로 lexical scopingenvironments를 사용할 때

model:

if logical: keep a copy of the model frame in the result? If the input value for model is a model frame (likely from an earlier call to the rpart function), then this frame is used rather than constructing new data.

공식이 까다로울 수있다. 테리 테 루노 (The Terry Therneau)는 model=TRUE과 함께이 두 패키지에 우리의 삶을 편하게 만들어 주셨습니다.

관련 문제