2017-05-19 24 views
3

成功/失敗データ(特定の期間に生き残り/死亡したツリー)があり、私の観測(7サイト)ごとに二項分布の誤差を推定したいと考えています。 。これまで私はglmを使っていました。成功/失敗のエラー推定R

s <- c(1,20,0,40,2,1,0) # success 
f <- c(2,0,20,4,50,0,1) # failure 

#for each observation I would calculate this error: 

error <- vector() 
z_scores <- vector() 
p_value <- vector() 

    for (i in 1:7) { 
    models <- glm (cbind (s[i], f[i]) ~ 1, family = 'binomial') 
    error [i] <- summary (models)$coefficients[2] 
    z_scores [i] <- summary (models)$coefficients[3] 
    p_value [i] <- summary (models)$coefficients[4] 
    } 

これは最善のアプローチでしょうか?

ここで二項分布の確率はどのように推定されますか?かかわらずsまたはfいずれかがここで=0

答えて

5

あるとき私の誤差が非常に高く、成功と失敗の数は使用せず(ゼロによって引き起こされる極端のものを除く)検索結果のほとんどを再計算するためにいくつかのコードです

注意glmの意味を説明し、その意味を説明します。

s <- c(1, 20, 0, 40, 2, 1, 0) # success 
f <- c(2, 0, 20, 4, 50, 0, 1) # failure 

#for each observation I would calculate this error: 

error <- vector() 
z_scores <- vector() 
p_value <- vector() 

for (i in 1:7) { 
    models <- glm(cbind(s[i], f[i]) ~ 1, family = 'binomial') 
    error[i] <- summary(models)$coefficients[2] 
    z_scores[i] <- summary(models)$coefficients[3] 
    p_value[i] <- summary(models)$coefficients[4] 
} 

logit <- function(x){ 
    log(x/(1 - x)) 
} 

dlogit <- function(x){ 
    1/x/(1 - x) 
} 

p_hat <- s/(s + f) 
## sqrt(p_hat * (1 - p_hat)/(s + f)) 
## is the standard error of p_hat 
## error1 is the standard error of logit(p_hat) 
error1 <- dlogit(p_hat) * sqrt(p_hat * (1 - p_hat)/(s + f)) 
## divide the estimation by the standard error, you get z-score 
z_scores1 <- logit(p_hat)/error1 
p_value1 <- 2 * pnorm(-abs(z_scores1)) 

あなたが知る必要がある最初の事は、我々が最初にこのケースでは、いくつかのモデルを(持っている、統計でなど標準誤差、Zスコア、p値の根拠とされ、二項モデル:s|(s+f) ~ Binomial(s + f, p))と我々は、データがランダムに生成されているので、我々は推定値である良い方法を知りたい、ここに来て)私たちが持っているデータと

1に合わせて、それを使用したい)推計(この場合はp

2を取得します標準誤差、zスコアとp値を「推定のランダム性を測定する」といいます。ここで重要な「トリック」があります: Tデータを生成し、真のメカニズムを知って、私たちは、およそa)は我々のモデルは、(またはそれに類似した)データの生成と

Bの真のメカニズムである仮定

を行うことにより、当社の推定におけるランダム性を計算することができます)実際のパラメータは私たちの見積もりに似ています(この場合、標本サイズはちょうどs + fであるので、s + fは推論(標準誤差、zスコアとp値)を検証するのに十分な大きさでなければなりません) 。そして、i = 1,6,7の場合、サンプルサイズは実際には小さく、対応する標準誤差、zスコア、p値が信じられないほどになることがわかります。

次に、私の計算の背後にある技術的な詳細とその意味について話すことができます。

logit(p) ~ N(mu, sigma^2)

をそしてロジット関数は、ちょうど私のコードではそのようなものです:glmで、Binomial(n, p)モデルのほかに、あなたもこのようなpするためのモデルを想定しています。この単純なケースで

は、二項確率 pの推定がちょうど p_hat <- s/(s + f)(使用 glmかどうか)で、二項変数の分散式から、我々は推定確率のための分散を得ることができます pここであれば、 p * (1 - p)/nですは、仮定bによって実際の pに似ており、 pを置き換えるために使用すると、推定誤差 pを得ることができます。CLTとデルタ法に従うと、サンプルサイズが十分に大きい場合、 s/(s + f)または logit(s/(s + f))を正規分布に従って扱うことができます。たとえば、 s/(s + f)は約 N(p, s * f/(s + f)^3)logit(s/(s + f))は約 N(logit(p), dlogit(s/(s + f))^2 * s * f/(s + f)^3)です。

単純に言えば、glmが計算する標準誤差、z-スコアおよびp-値は、logit(s/(s + f))の標準誤差、z-スコアおよびp-値です。これらは帰無仮説の有効な結果であるlogit(p) = 0、つまりp = 0.5です。したがって、glmから得られたz-スコアおよびp-値は、s + fが大きい場合にsおよびfが等しい確率で起こるかどうかをテストすることです。

そして、私は0 sまたはfは0、fまたはsが、これが本当であれば、1になりますが起こるの推定確率等しいによって引き起こされる極端な値について話しますが、データの発生機構は、実際には非ランダムであります!!当初、私たちは推定値を使って推定値をランダムに計算し、sまたはfが0に等しいと仮定して、推定値を真実として使用すると、100%これはちょっとばかげている。そのような場合、glmのような多くの方法は有効ではありません。一般に、サンプルサイズがs + fの場合、s = 0またはf = 0の場合、sまたはfの確率は実際には小さいと考えられますが、ケース6または7のようにサンプルサイズが実際には小さければ、結論。要するに

二項モデルはglm結果から、真であるならば、私のコードと私の分析は、上記のように、我々はケースi = 2, 3, 4, 5には、sfの確率が大幅に互いに異なっていると言うことができます。

関連する問題