So, da bin ich nochmal.
Aaalso: Bei der normalen OLS-Regression erhälst Du super schnell einen Vorschlag für den Intercept und einen für die Steigung, beide mit Standardfehler, sodass man daraus vielleicht einen Test bauen könnte - wie gesagt, kann vielleicht ein anderer beobachten.
Bei der Bayes-Regression geht alles viel langsamer, dafür erhälst Du dann aber nicht nur einen Wert sondern eine Stichprobe aus der Verteilung möglicher Regressionsgeraden. Im Beispiel oben also anstelle jeder linearen Regressionsgeraden einen Dataframe der 4000 mögliche Regressiongeraden beschreibt. Natürlich kannst Du mehr als 4000 haben, wenn Du mehr Rechenzeit inkauf nimmst. Mit diesen Regressionsgeraden lässt sich sehr viel auswerten.
Um das für jedes Land getrennt zu rechnen habe ich das in eine for-Schleife gesetzt.
Ich habe mal eine Funktion compare geschrieben, mit der ich an einem einzelnen Land vormache, wie man die Werte für den Intercept und für den ersten Koeffizienten extrahiert und in Grafiken darstellt:
Code: Alles auswählen
compare <- function(land){
par(mfrow = c(2, 2))
i1 <- as.data.frame(out1[[land]])$`(Intercept)`
i2 <- as.data.frame(out2[[land]])$`(Intercept)`
y1 <- as.data.frame(out1[[land]])$y1
y2 <- as.data.frame(out2[[land]])$y2
mii <- min(c(i1, i2))
mai <- max(c(i1, i2))
miy <- min(c(y1, y2))
may <- max(c(y1, y2))
hist(i1, breaks = 30, xlim = c(mii, mai), main = "Intercept to y1",
xlab = "Intercept")
hist(y1, breaks = 30, xlim = c(miy, may), main = "Coefficent for y1",
xlab = expression(beta[1]))
hist(i2, breaks = 30, xlim = c(mii, mai), main = "Intercept to y2",
xlab = "Intercept")
hist(y2, breaks = 30, xlim = c(miy, may), main = "Coefficent for y2",
xlab = expression(beta[1]))
cat("\ninvestigating land ");cat(land);cat(":\n")
cat("intercept 1 > intercept 2 in ");cat(sum(i1 > i2)/4000 * 100);cat(" %\n")
cat("y1 > y2 in ");cat(sum(y1 > y2)/4000 * 100);cat(" %\n")
}
Wenn wir das jetzt mal für Land 4 machen:
dann erhalten wir vier Histogramme und zwei zahlenmäßige Vergleiche.
In der linken Spalte sehen wir, dass der mögliche Intercept für Land 4 im einen Fall viel weniger streut als im anderen. Dafür ist in einem Fall unklar, ob der Koeffizient für y1 positiv oder negativ sein soll mit Tendenz zum negativen, im anderen ist eher klar, dass der Koeffizient für y2 dazu tendiert, größer als Null zu sein:
- Rplot02.png (9.07 KiB) 640 mal betrachtet
Insgesamt überschneiden sich die Koeffizienten in beiden Fällen aber deutlich. Die Berechnung sagt uns, dass der Koeffizient von y1 nur in 7,5% größer als der von y2 ist. In der Bayes-Statistik kann man sagen: Der y2 ist mit 92,5% Wahrscheinlichkeit größer als der Koeffizient für y1. In der klassischen NHST-Statistik wäre so ein Satz nicht möglich.
Probier mal compare mit den anderen Zahlen zwischen 1 und 6 aus. Wirklich interessant wird es natürlich erst mit den richtigen Daten.
So, ich habe keine Ahnung, ob das eine Richtung ist, in die Du gehen kannst oder willst, drum höre ich jetzt mit dieser Richtungsangabe und den Stichworten "Bayes-Statistik" und "rstanarm" mal auf. Wenn Du willst, kannst Du gerne nachfragen. Sonst hoffe ich für Dich, dass noch jemand erklärt, wie man das mit Mitteln der klassichen Regressionrechnung sauber macht.
LG,
Bernhard