2

バイナリロジットモデルで「もしも」の状況をモデル化しようとしています。テストの難易度 (1 = 最も簡単、5 = 最も難しい) を考慮して、性別を対照として、テストに合格する確率を推定しています。（データはこちら）。学生は、一般的に厳しいテストを受けます (データの「高」)。これから、合格の可能性に対するテストの難易度の影響を推定できます。

model = glm(PASS ~ as.factor(SEX) + as.factor(HIGH), family=binomial(link="logit"), data=df)
summary(model)

また、次の方法で合格の予測確率を取得することもできます。

predict.high = predict(model, type="response")

問題は、代わりに「LOW」テストが与えられた場合はどうなるかということです。新しい確率を取得するには、次のようにします。

newdata = rename.vars(subset(df, select=c(-HIGH)), 'LOW','HIGH')
predict.low = predict(model, newdata=newdata, type="response")

しかし、この場合、何人の追加の学生が合格したかをどのように知ることができますか? 私が見ていないのに明らかなスイッチはありglm()ますか？

4

2 に答える 2

3

私は、Gelman と Hill (2006) に基づいて書いた予測用のコードをまだ掘り下げようとはしていません。今でもそうするつもりです。私の限られた経験の中であなたの質問のユニークな側面の 1 つは、私が 1 つの観察 (この場合は 1 つのテストを受ける 1 人の学生) を予測することに慣れていたことです。しかし、あなたは 2 つの予測セットの違いを予測したいと考えているようです。言い換えると、5 つの難しい試験ではなく 5 つの簡単な試験を与えられた場合に合格する学生の数を予測したいと考えています。

Gelman と Hill (2006) がそれをカバーしているかどうかはわかりません。また、頻度主義的なアプローチでこれを行いたいようです。

単一の観測値を予測して、各観測値の信頼区間を得ることができれば、おそらく各グループ内の加重平均通過確率を推定し、2 つの加重平均を差し引くことができると考えています。デルタ法を使用して、加重平均とその差の信頼区間を推定できます。

そのアプローチを実装するには、予測された観測値間の共分散を 0 と仮定する必要がある場合があります。

共分散が 0 であると仮定しても満足できない場合は、おそらくベイジアン アプローチの方が適しています。繰り返しますが、私は単一の観測値の予測にしか慣れていません。ベイジアン アプローチを使用して、独立変数を含めて 1 つの観測値を予測しましたが、従属変数は含めずに、観測値を予測しました。同じベイジアン実行ですべての観測値を予測できると思います (各学生を HIGH と LOW で予測します)。各グループの合格テストの加重平均と加重平均の差は派生パラメーターであり、ベイズロジスティック回帰のコードに直接含めることができると思います。次に、テストの各グループに合格する確率と、テストの各グループに合格する確率の差について、ポイント推定値と分散推定値を取得します。

これまでのところ、この回答は必要以上に会話的なものでした。私は、それらの戦略を実行しようとする時間がないまま、試みるべき戦略を単に計画しているだけです。R と WinBUGS のすべてのコードを提供して、提案された両方の戦略を実装するには、数日かかるかもしれません。(WinBUGS または OpenBUGS は R 内から呼び出すことができます。) 進行するにつれて、この回答にコードを追加します。私の提案した戦略や今後のコードが間違っていると思われる人がいたら、遠慮なく間違いを指摘して訂正してくれることを願っています。

編集

以下は、フェイク データを生成し、頻度論的アプローチとベイジアン アプローチを使用してそのデータを分析するコードです。上記の予測のアイデアを実装するためのコードはまだ追加していません。今後 1 ～ 2 日でベイジアン予測コードを追加しようと思います。5 つではなく 3 つのテストのみを使用しました。以下のコードの書き方では、生徒数 n を 6 つの等しい整数に分割できるゼロ以外の任意の数に変更できます。

# Bayesian_logistic_regression_June2012.r
# June 24, 2012

library(R2WinBUGS)
library(arm)
library(BRugs)

set.seed(3234)


# create fake data for n students and three tests

n <- 1200

# create factors for n/6 students in each of 6 categories

gender <- c(rep(0, (n/2)), rep(1, (n/2)))
test2  <- c(rep(0, (n/6)), rep(1, (n/6)), rep(0, (n/6)),
            rep(0, (n/6)), rep(1, (n/6)), rep(0, (n/6)))
test3  <- c(rep(0, (n/6)), rep(0, (n/6)), rep(1, (n/6)),
            rep(0, (n/6)), rep(0, (n/6)), rep(1, (n/6)))

# assign slopes to factors

B0      <-  0.4
Bgender <- -0.2
Btest2  <-  0.6
Btest3  <-  1.2

# estimate probability of passing test

p.pass <- (     exp(B0 + Bgender * gender + 
                         Btest2  * test2  + 
                         Btest3  * test3) /
           (1 + exp(B0 + Bgender * gender +
                         Btest2  * test2  + 
                         Btest3  * test3)))

# identify which students passed their test, 0 = fail, 1 = pass

passed   <- rep(0, n)
r.passed <- runif(n,0,1)
passed[r.passed <= p.pass] = 1

# use frequentist approach in R to estimate probability
# of passing test

m.freq <- glm(passed ~ as.factor(gender) +
                       as.factor(test2)  +
                       as.factor(test3)  , 
                       family = binomial)
summary(m.freq)

# predict(m.freq, type = "response")


# use OpenBUGS to analyze same data set

# Define model

sink("Bayesian.logistic.regression.txt")
cat("
model {

# Priors

 alpha ~ dnorm(0,0.01)
 bgender ~ dnorm(0,0.01)
 btest2 ~ dnorm(0,0.01)
 btest3 ~ dnorm(0,0.01)

# Likelihood

 for (i in 1:n) {
    passed[i] ~ dbin(p[i], 1)
    logit(p[i]) <- (alpha + bgender * gender[i] +
                            btest2  * test2[i]  +
                            btest3  * test3[i])
 }

# Derived parameters

 p.g.t1 <- exp(alpha) / (1 + exp(alpha))
 p.b.t1 <- exp(alpha + bgender) / (1 + exp(alpha + bgender))

 p.g.t2 <- (    exp(alpha +           btest2) / 
           (1 + exp(alpha +           btest2)))
 p.b.t2 <- (    exp(alpha + bgender + btest2) / 
           (1 + exp(alpha + bgender + btest2)))

 p.g.t3 <- (    exp(alpha +           btest3) / 
           (1 + exp(alpha +           btest3)))
 p.b.t3 <- (    exp(alpha + bgender + btest3) / 
           (1 + exp(alpha + bgender + btest3)))

}

", fill = TRUE)
sink()

my.data <- list(passed = passed, 
                gender = gender,
                test2  = test2,
                test3  = test3, 
                n      = length(passed))

# Inits function

inits <- function(){ list(alpha   = rlnorm(1), 
                          bgender = rlnorm(1),
                          btest2  = rlnorm(1),
                          btest3  = rlnorm(1)) }

# Parameters to estimate

params <- c("alpha", "bgender", "btest2", "btest3", 
            "p.g.t1", "p.b.t1", "p.g.t2", "p.b.t2",
            "p.g.t3", "p.b.t3")

# MCMC settings

nc <- 3
ni <- 2000
nb <- 500
nt <- 2

# Start Gibbs sampling

out <- bugs(data = my.data, inits = inits,
parameters.to.save = params, 
"c:/users/Mark W Miller/documents/Bayesian.logistic.regression.txt",
program = 'OpenBUGS', 
n.thin = nt, n.chains = nc, 
n.burnin = nb, n.iter = ni, debug = TRUE)

print(out, dig = 5)

予測に加重平均アプローチを実装しようとする前に、それが機能する可能性があることを確信したかったのです。そこで、次のコードを作成しました。

# specify number of girls taking each test and
# number of boys taking each test

g.t1 <- rep(0,400)
b.t1 <- rep(0,120)
g.t2 <- rep(0,1200)
b.t2 <- rep(0,50)
g.t3 <- rep(0,1000)
b.t3 <- rep(0,2000)

# specify probability of individuals in each of the
# 6 groups passing their test

p.g1.t1 <- 0.40
p.b1.t1 <- 0.30
p.g1.t2 <- 0.60
p.b1.t2 <- 0.50
p.g1.t3 <- 0.80
p.b1.t3 <- 0.70

# identify which individuals in each group passed their test

g.t1[1:(p.g1.t1 * length(g.t1))] = 1
sum(g.t1)

b.t1[1:(p.b1.t1 * length(b.t1))] = 1
sum(b.t1)

g.t2[1:(p.g1.t2 * length(g.t2))] = 1
sum(g.t2)

b.t2[1:(p.b1.t2 * length(b.t2))] = 1
sum(b.t2)

g.t3[1:(p.g1.t3 * length(g.t3))] = 1
sum(g.t3)

b.t3[1:(p.b1.t3 * length(b.t3))] = 1
sum(b.t3)

# determine the weighted average probability of passing
# on test day for all individuals as a class

wt.ave.p <- ((p.g1.t1 * length(g.t1) + p.b1.t1 * length(b.t1) +
 p.g1.t2 * length(g.t2) + p.b1.t2 * length(b.t2) +
 p.g1.t3 * length(g.t3) + p.b1.t3 * length(b.t3) ) / 

 (length(g.t1) + length(b.t1) + length(g.t2) + 
  length(b.t2) + length(g.t3) + length(b.t3)))

wt.ave.p

# determine the expected number of individuals passing
# their test in the class as a whole

exp.num.pass <- wt.ave.p *  (length(g.t1) + length(b.t1) +
                             length(g.t2) + length(b.t2) +
                             length(g.t3) + length(b.t3))
exp.num.pass

# determine the number of individuals passing

num.passing <- (sum(g.t1) + sum(b.t1) + 
                sum(g.t2) + sum(b.t2) + 
                sum(g.t3) + sum(b.t3) )
num.passing

# the expected number of students passing, exp.num.pass,
# should equal the observed number of students passing,
# num.passing regardless of the number of students in each
# group and regardless of the probability of passing a 
# given test, within rounding error

identical(round(exp.num.pass), round(num.passing))

うまくいけば、次の数日で、予測コードを上記のベイジアン コードに追加してみることができます。

編集 - 2012 年 6 月 27 日

私はこれを忘れていません。むしろ、いくつかの問題に遭遇しました。

  1. ロジスティック回帰を使用すると、a) 特定のグループの学生がテストに合格する確率 p、および b) 特定の学生がテストを受ける結果 (0 または 1) を予測できます。次に、すべての 0 と 1 が平均化されます。これらのどれを使用するかわかりません。予測された p の点推定値と SD は、既知のテスト結果の推定 p と同じです。予測された 0 と 1 の平均の点推定値は少し異なり、平均化された 0 と 1 の SD ははるかに大きくなります。予測された 0 と 1 の平均である b が必要だと思います。ただ、念のため、いろいろなサイトや本を調べてみています。コレット (1991) は、コンピューター コードを使用しない実例を示しています。

  2. 多くの派生パラメーターを使用すると、プログラムの実行に時間がかかります。

  3. どうやら OpenBUGS は、予測コードがなくても頻繁にクラッシュしているようです。それが私が間違ったことをしているためなのか、Rの最近のバージョンでの変更やRパッケージの最近のバージョンでの変更によるものなのか、それとも64ビットのRなどでコードを実行しようとしているからなのかはわかりません。そうしないと。

すぐに予測コードを投稿しようとしますが、上記の問題のすべてが私の速度を低下させました。

于 2012-06-24T14:58:58.617 に答える
0

このアプローチを使用して、カットオフを簡単に見つけることができます。

cutoff <- runif(length(predicted_probabilities))

これは、Metropolis-Hastings に基づく決定論的な決定です。

于 2014-10-08T21:03:05.820 に答える