r - 如何从 R 中的包中提取信息并在函数中使用？

Question

对于模糊的问题标题，我深表歉意。我想做的是使用 geepack R 包中的 geeglm 在 R 中运行回归，然后使用其中的信息来计算准似然信息标准(QIC；Pan 2001)。我可以很容易地为单个模型做到这一点，但我想编写一个通用函数来为各种不同类型的模型做到这一点。我想我真正的问题是是否有比一长串嵌套的 ifelse 语句更好的选择？

这是我当前的代码:

library(geepack)
data(dietox) #data from the geepack package
# Run gee regression
dietox$Cu <- as.factor(dietox$Cu)
mf <- formula(Weight ~ Cu * (Time + I(Time^2) + I(Time^3)))
gee1 <- geeglm(mf, data = dietox, id = Pig, family = gaussian, corstr = "ar1")

然后我可以运行一个函数来计算拟似然性:

QlogLik.normal <- function(model.R) {
  library(MASS)
  mu.R <- model.R$fitted.values
  y <- model.R$y
  # Quasi Likelihood for Normal
  quasi.R <- sum(((y - mu.R)^2)/-2)
  quasi.R
  }

但是，我想编写一个更通用的函数，因为拟似然函数对于每个分布都是不同的。上述函数适用于 gee1，因为它具有高斯(正态)分布。如果我想将它概括为各种发行版，我可以使用一系列嵌套的 ifelse 语句(如下)，但我不知道这是否是最好的方法。有没有人有其他选择或更好的解决方案？这至少可以说不是很优雅(显然我没有太多编程或 R 经验)。

QlogLik <- function(model.R) {
  library(MASS)
  mu.R <- model.R$fitted.values
  y <- model.R$y
  ifelse(model.R$modelInfo$variance == "poisson",
     # Quasi Likelihood for Poisson
     quasi.R <- sum((y*log(mu.R)) - mu.R),
     ifelse(model.R$modelInfo$variance == "gaussian",
       # Quasi Likelihood for Normal
       quasi.R <- sum(((y - mu.R)^2)/-2),
       ifelse(model.R$modelInfo$variance == "binomial",
         # Quasilikelihood for Binomial
         quasi.R <- sum(y*log(mu.R/(1 - mu.R)) + log(1 - mu.R)),
         quasi.R <- "Error: distribution not recognized")))
  quasi.R
  }

在此示例中，我使用 geeglm 的模型输出来提取用于对方差建模的分布类型

 model.R$modelInfo$variance

但可能还有其他方法可以确定 geeglm 模型中使用的分布。任何帮助将不胜感激。

Answer 1

您应该能够像这样重写您的函数:

QlogLik <- function(model.R) {
  library(MASS)
  mu.R <- model.R$fitted.values
  y <- model.R$y
  type <- family(model.R)$family
  switch(type,
         poisson = sum((y*log(mu.R)) - mu.R),
         gaussian = sum(((y - mu.R)^2)/-2),
         binomial = sum(y*log(mu.R/(1 - mu.R)) + log(1 - mu.R)),
         stop("Error: distribution not recognized"))
}

正如@baptise 指出的那样，switch 在这些情况下很有用。您使用 family(model.R)$family 来自动检测应该与 switch 一起使用的族类型。

此外，如果您在不同情况下执行的命令超过一行，您可以用大括号将这些行括起来({ do something here })。

switch(type,
       type1 = { something <- do(this)
                 thisis(something) },
       type2 = do(that))                      

希望对您有所帮助!