函數內部公式錯誤

Question

我想在函數內部使用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) 

Answer 1

這是作用域的問題。 注意的basehaz環境：

environment(basehaz) 
<environment: namespace:survival> 

同時：

environment(BaseFun) 
<environment: R_GlobalEnv> 

所以這就是爲什麼功能basehaz找不到函數內部的局部變量。

一個可能的解決方案是發送X到使用assign頂部：

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) 

其他解決方案可能涉及的環境更直接打交道。

Answer 2

我正在追蹤@ moli對@ aatrujillob的回答的評論。他們是有幫助的，所以我想我會解釋它是如何爲我解決的，以及與rpart和partykit包類似的問題。

一些玩具數據：

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對象不存在。

my_formula <- Trt~W+X 
exists("train_data") 
# [1] FALSE 
exists("train_data", envir = environment(my_formula)) 
# [1] FALSE 

這是我的功能，它與原始功能類似。再次，

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()就像是coxph()中，你可以指定參數model=TRUE來解決這些問題。這裏再次提出這個額外的論點。

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 

文檔model論點rpart()：因爲他們的方式，並不總是天然的（對我來說）使用lexical scoping和environments

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.

公式可能會非常棘手。謝天謝地Terry Therneau讓我們的生活變得更輕鬆，model=TRUE在這兩個包裏！