これは私の最初の投稿なので、改善すべき点や追加情報があればお知らせください。
10.11.6 を実行している Mac で R 3.4.3 を使用しています。
カウントデータを扱っています。負の二項分布のシータ/分散/サイズ パラメーター (k) に興味があります (これらの用語は同じ意味で使用されていると理解しています)。fitdistrplus
関数を使用して、パッケージを使用して最尤法でパラメーターを推定するこれらのデータに NB 分布を当てはめていfitdist
ます。データのさまざまなグループが、さまざまな分布を使用してより適切にモデル化されているかどうかに興味があります。したがって、分布をすべてのデータに当てはめます。次に、データは 1 つの 2 水準因子に基づいて分割され、分布はこれら 2 つの別々のグループに当てはまります。
分布をデータ セット全体に当てはめると、標準誤差を含む mu とサイズの推定値が得られます。次に、データを分割します。この同じアプローチは、データの半分 (グループ A) ではうまく機能しますが、理論的には構造的に同じであるはずの残りの半分 (グループ B) では機能しません。代わりに、mu とサイズの推定値を取得しますが、標準誤差の NA を使用します。
optim
背後の関数にlower = c(0,0) と upper = c(inf, inf) を課してfitdist
も何も達成されず、出力が NaN ではなく NA であるため、しようとしているのとは何の関係もないと思いますとにかく負の数を見積もってください(よく議論されるエラー100)。
そして、それが何か関係があり、それが何もしなかった場合に備えて、私はすべてのゼロを削除しました。
私の質問は、なぜ NA が生成されるのか (そして、最終的に推定値の標準誤差を取得するにはどうすればよいのか) です。
ここに私のデータと私のコードがあります:
require(fitdistrplus)
data.set <- read.csv("data.set.csv")
count.A <- subset(data.set, category == "A")
count.B <- subset(data.set, category == "B")
# All Data
plotdist(data.set$count, histo = TRUE, demp = TRUE)
count.nb <- fitdist(data.set$count, "nbinom")
plot(count.nb)
LL.nb <- logLik(count.nb)
count.p <- fitdist(data.set$count, "pois")
plot(count.p)
LL.p <- logLik(count.p)
cdfcomp(list(count.p, count.nb),legendtext = c("Poisson", "negative binomial"))
gofstat(list(count.p, count.nb),fitnames = c("Poisson", "negative binomial"))
# Group A
plotdist(count.A$count, histo = TRUE, demp = TRUE)
A.count.nb <- fitdist(count.A$count, "nbinom")
plot(A.count.nb)
A.LL.nb <- logLik(A.count.nb)
A.count.p <- fitdist(count.A$count, "pois")
plot(A.count.p)
A.LL.p <- logLik(A.count.p)
cdfcomp(list(A.count.p, A.count.nb),legendtext = c("Poisson", "negative binomial"))
gofstat(list(A.count.p, A.count.nb),fitnames = c("Poisson", "negative binomial"))
# Group B
plotdist(count.B$count, histo = TRUE, demp = TRUE)
B.count.nb <- fitdist(count.B$count, "nbinom", method = "mle")
plot(B.count.nb)
B.LL.nb <- logLik(B.count.nb)
B.count.p <- fitdist(count.B$count, "pois")
plot(B.count.p)
B.LL.p <- logLik(B.count.p)
cdfcomp(list(B.count.p, B.count.nb),legendtext = c("Poisson", "negative binomial"))
gofstat(list(B.count.p, B.count.nb),fitnames = c("Poisson", "negative binomial"))