与生存包的 R 的 logsum 函数(mlogit 包)类似的功能

Similar function to R's logsum function (mlogit package) for Survival package

我试图使用 R 的 survival 包(clogit 函数)为每个人获得预期的效用,但我无法找到一个简单的解决方案,例如 mlogit 的 logsum .

下面我举了一个例子,说明如何使用 mlogit 包。它非常简单:它只需要使用 mlogit 函数对变量进行回归,保存输出并将其用作 logsum 函数中的参数——如果需要,这里有一个简短的解释vignette. And what I want is to know the similar method for clogit. I've read the package's manual 但我没能掌握执行分析的最适当功能。

注意 1:我偏爱像 mlogit's 这样的函数是因为我以后可能需要执行大量的回归,并且能够在不同的场景中执行正确的估计会很有帮助。

注意 2:我不希望下面创建的数据集代表数据的行为方式。我设置这个例子只是为了在逻辑回归之后执行函数。

**

library(survival)
library(mlogit)

#creating a dataset

df_test=data.frame(id=rep(1:20,each=4),
                   choice=rep(c("train","car","plane","boat")),
                   distance=c(rnorm(80)*10),
                   )

f=function(x,y,z) {
    
  v=round(rnorm(x,y,z))
    
    while(sum(v)>1 | sum(v)==0) {
      
      v=round(rnorm(x,y,z))
      
    }
  
return(v)
    
}

result1=c()

for (i in 1:20) {
  
  result=f(4,0.5,0.1)
  
  result1=c(result,result1)
  
}

df_test$distance=ifelse(df_test$distance<0,df_test$distance*-1,df_test$distance)
df_test$price = 0
df_test$price[df_test$choice=="plane"] = rnorm(20, mean = 300, sd=30)
df_test$price[df_test$choice=="car"] = rnorm(20, mean = 50, sd=10)
df_test$price[df_test$choice=="boat"] = rnorm(20, mean = 100, sd=15)
df_test$price[df_test$choice=="train"] = rnorm(20, mean = 120, sd=25)

df_test$choice2=result1
           
mlog=mlogit(choice2 ~ distance + price , data = df_test)

#the function logsum generates expected utility for each individual

logsum(mlog)

#so what would be adequate alternative with survival's clogit? I set an exemple below of
#of what i would like to regress and then perform something like logsum()

clog=clogit(choice2 ~ distance + price + as.factor(choice), strata(id), data = df_test)

**

你提供的小插图说对数计算如下:

根据我的阅读,这类似于用于构建“线性预测器”的计算。 lp 是 t(coef(clog)) %*% Xhat。如果我的解释是正确的,那么它存储在 clog-object:

 clog["linear.predictor"]

所以你可以采取:

log(sum( exp(clog[["linear.predictors"]]) ))
[1] 4.286592

如果您需要按 df_test$choice 进行分隔,则(因为 clog 的 'linear.predictors' 元素的长度与 'choice' 列的长度相同df_test) 只是:`

 tapply(clog$linear.predictors, df_test$choice, function(x){ 
               log(sum( exp(x) ))})
 #--------------------------------------
    boat      car    plane    train 
2.824502 3.506756 2.825004 1.734258 

如果您需要按 id 进行汇总,您可以这样做:

tapply(clog$linear.predictors, df_test$id, function(x){ 
  log(sum( exp(x) ))})

       1        2        3        4        5        6        7        8        9 
1.405896 1.506152 1.454507 1.539060 1.467082 1.428482 1.393582 1.521877 1.480670 
      10       11       12       13       14       15       16       17       18 
1.466490 1.416338 1.500581 1.528075 1.488253 1.398405 1.445014 1.483623 1.404041 
      19       20 
1.460672 1.452709