2

我为模糊的问题标题道歉。我想做的是使用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 模型中使用的分布。任何帮助，将不胜感激。

4

2 回答 2

4

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

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))

我希望这有帮助！

于 2012-09-16T08:13:35.613 回答
2

您也可以使用model.R$family$familywhich 给出用于模拟方差的分布类型，但到目前为止我不知道您是否可以消除这些ifelse陈述。您的quasi.R代码中的 不同发行版之间有所不同，因此您必须分别定义它们中的每一个。

顺便说一句，这是一个很好的问题，感谢您发布它：我过去也遇到过类似的情况，希望得到一些关于如何更有效地编写代码的建议。

于 2012-09-16T04:01:56.583 回答