R: Vorhersage (0,1) in die Logistische regression in glm()

Habe ich noch nicht versucht zu Graben, mein code für die Vorhersage, dass ich schrieb, basierend auf Gelman und Hill (2006), die, ich glaube mich zu erinnern simulation. Ich habe immer noch die Absicht, das zu tun. Ein Aspekt Ihrer Frage, schien einzigartig in meiner begrenzten Erfahrung war, die ich gewohnt war, um die Vorhersage für eine einzelne Beobachtung (in diesem Fall ein einziger Schüler ein Einzel-test). Sie scheinen jedoch Vorhersagen möchten, einen Unterschied zwischen zwei Datensätzen Vorhersagen. In anderen Worten, Sie wollen Vorhersagen, wie viele weitere Schüler übergeben, wenn gegeben, eine Reihe von 5 einfachen Prüfungen eher als die ein-Satz von 5 harten Prüfungen.

Ich bin mir nicht sicher, ob Gelman und Hill (2006) abgedeckt. Sie auch scheinen möchten, tun Sie dies mit einem frequentistischen Ansatz.

Ich denke, dass, wenn Sie Vorhersagen kann für eine einzelne Beobachtung, so dass Sie ein Konfidenzintervall für jede Beobachtung, dann kann man vielleicht schätzen die gewichtete Durchschnittliche Wahrscheinlichkeit der Weitergabe innerhalb jeder Gruppe und subtrahieren Sie die beiden gewichtete Durchschnitte. Die delta-Methode verwendet werden könnte, um schätzen eines Konfidenzintervalls auf der gewichteten Mittelwerte und deren Differenz.

Kovarianz zwischen vorhergesagten Beobachtungen kann davon ausgegangen werden, auf 0 zu implementieren, der Ansatz.

Angenommen, wenn eine Kovarianz von 0 ist nicht zufriedenstellend, dann vielleicht ein Bayes-Ansatz besser wäre. Wieder, ich bin nur vertraut mit der Vorhersage für eine einzelne Beobachtung. Mit einem Bayes-Ansatz habe ich vorausgesagt eine einzige Beobachtung, indem die unabhängigen Variablen, nicht aber die abhängige variable, für die Beobachtung zu prognostizieren. Ich nehme an, Sie könnten Vorhersagen für jede Beobachtung in der gleichen Bayes-run (Vorhersagen, jeder Schüler in der HIGH und LOW). Die gewichteten Durchschnitte der übergabe von tests für jede Gruppe und die Differenz in den gewichteten Durchschnitt sind die abgeleiteten Parameter und ich vermute, dass könnte sich direkt in den code für die Bayes-logistic regression. Dann würden Sie Ihr Punkt-Schätzung und die Schätzung der Varianz für die Wahrscheinlichkeit der Weitergabe jede Gruppe von tests für den Unterschied in der Wahrscheinlichkeit der Weitergabe jede Gruppe von tests. Wenn Sie möchten, dass der Unterschied in der Anzahl der Studenten bestehen jede Gruppe von tests, die vielleicht aufgenommen werden könnte in die Bayes-code als abgeleitete parameter auch.

Erkenne ich diese beantworten, so weit, hat mehr Konversation als gewünscht werden könnten. Ich bin einfach der mapping-Strategien zu versuchen, ohne die Zeit hatte, noch zu versuchen, die Umsetzung dieser Strategien. Alle R und WinBUGS-code zu implementieren, der sowohl die vorgeschlagenen Strategien kostet mich vielleicht ein paar Tage. (WinBUGS oder OpenBUGS aufgerufen werden kann, innerhalb R.) ich werde fügen Sie den code, um diese Antwort, als ich entlang gehen. Wenn jemand hält meine vorgeschlagenen Strategien und/oder bevorstehende code, falsch, ich hoffe, Sie fühlen sich frei, um darauf hin meine Fehler und Korrekturen bieten.

BEARBEITEN

Unten ist code, der generiert gefälschte Daten und analysiert diese Daten mit einem frequentistischen und Bayes-Ansatz. Ich habe nicht noch den code zum implementieren der oben genannten Ideen für die Vorhersage. Ich werde versuchen, fügen Sie die Bayes ' sche Vorhersage code in den nächsten 1-2 Tagen. Ich habe nur drei tests statt fünf. Die Art und Weise Sie den code unten geschrieben Sie können ändern, die Anzahl der Studierenden, n, um alle nicht-null-Zahl, die geteilt werden können in 6 gleichen ganzen zahlen.

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

Bevor ich versucht zu implementieren, der gewichtete-Durchschnitt-Ansatz zur Vorhersage ich wollte mich selbst davon zu überzeugen, dass es funktionieren könnte. Also ich ginned bis der folgende code, der zu zeigen scheint, kann Sie:

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

Hoffentlich in den nächsten paar Tagen kann ich versuchen die Vorhersage-code zu den oben genannten Bayes-code.

BEARBEITEN - Juni 27, 2012

Habe ich auch nicht vergessen zu diesem. Ich habe stattdessen traf auf mehrere Probleme:

1. Logistischen regression ist es möglich, vorherzusagen: a) die Wahrscheinlichkeit p, dass die Schüler in einer bestimmten Gruppe einen test, und b) dem Ergebnis eines bestimmten Schüler einen test (0 oder 1). Alle die 0 und 1 werden dann gemittelt. Ich bin mir nicht sicher, welche dieser zu nutzen. Die Punkt-Schätzung und SD der vorhergesagten p ist identisch mit der geschätzten p nach bekannten Untersuchungsergebnissen. Die Punkt-Schätzung der durchschnittlichen vorhergesagten 0 s und 1 s ist ein wenig anders und der SD der Durchschnitt 0 ist, und 1 ist viel größer. Ich glaube, ich will b, der Mittelwert der vorhergesagten 0 und 1 ist. Allerdings bin ich versucht, zu untersuchen, auf verschiedenen websites und Bücher, um sicher zu sein. Collett (1991) ist ein Beispiel, dass nicht beschäftigen computer-code, aber das klappte Beispiel umfasst ein halbes Dutzend Variablen, einschließlich 2 Interaktionen, und ich bin mit ein wenig Mühe, mein Bayes-Schätzungen, passend zu Ihren frequentistischen Schätzungen.

2. Mit vielen abgeleiteten Parameter des Programms ist es, eine lange Zeit zu laufen.

3. Offenbar OpenBUGS wurde Häufig abstürzt, glaube ich, auch ohne Vorhersage-code. Ich bin mir nicht sicher, ob das ist, weil das, was ich falsch mache oder aufgrund von änderungen in den neueren Versionen von R oder änderungen in neueren Versionen von R-Pakete oder vielleicht auch, weil ich versuche, um den code auszuführen, der mit einem 64-bit-R oder etwas anderes.

Werde ich versuchen, nach der Vorhersage von code schnell, aber alle der oben genannten Fragen haben sich verlangsamt mich nach unten.