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ます。

また、さまざまなケースで実行するコマンドが 1 行を超える場合は、代わりに中括弧 ( { do something here }) で行を折り返すことができます。

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

これが役立つことを願っています！

于 2012-09-16T08:13:35.613 に答える
2

分散をモデル化するために使用される分布のタイプを与える which を使用することもできますが、これまでのところ、それらのステートメントmodel.R$family$familyを削除できるかどうかはわかりませんでした。ifelseコード内のquasi.Rはディストリビューションによって異なるため、それぞれを個別に定義する必要があります。

ところで、それは良い質問であり、投稿してくれてありがとう: 私は過去に同様の状況にあったので、コードをより効率的に記述する方法についてアドバイスを得たいと思っています。

于 2012-09-16T04:01:56.583 に答える