Gleichheit von Koeffizienten ermitteln

Wie rufe ich R-Funktionen auf, wie selektiere ich Daten, ich weiß nicht genau ....

Moderatoren: EDi, jogo

Antworten
h.itup
Beiträge: 11
Registriert: So Jun 06, 2021 5:34 pm

Gleichheit von Koeffizienten ermitteln

Beitrag von h.itup »

Guten Abend,

für einen Datensatz führe ich zwei verschiedene Regressionen für jedes Land durch.

Code: Alles auswählen

output1 <- lme4::lmList(x1 ~ y1| land, data = datensatz) 
output2 <- lme4::lmList(x2 ~ y2 | land, data = datensatz) 
Anschließend möchte ich ermitteln ob die Regressionskoeffizienten von Output 2 sich signifikant von den durchschnittlichen Koeffizienten der anderen Regression aus Output 2 auf Basis der Gruppenzugehörigkeit unterscheiden. Ich möchte bspw. wissen, ob für Gruppe 20 der Reg.Koeffizient aus Output 1 gleich zum Reg.Koeffizient der Gruppe 20 aus Output 2.

Ich scheitere leider daran mir einen solchen Befehl zu entwickeln. Ich habe einen Beispieldatensatz beigefügt. Der echte Datensatz ist noch einmal größer und beinhaltet mehr Länder sowie Gruppen.

Ich hoffe ich bin hier im Forum richtig.
Danke im Voraus.
Dateianhänge
datensatz.csv
(1.23 KiB) 77-mal heruntergeladen
bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: Gleichheit von Koeffizienten ermitteln

Beitrag von bigben »

Hallo!
h.itup hat geschrieben: Mi Jul 21, 2021 5:27 pmIch hoffe ich bin hier im Forum richtig.
Im Prinzip ja, aber der Sommer ist heiß, es sind Ferien, hier ist gerade wenig los. Ich kann Dir damit direkt leider nicht helfen und es kann sein, dass es ein paar Tage dauert, bis ein kompetenterer kommt.

Ich persönlich würde das Problem so angehen, dass ich anstelle der OLS-Regression eine Bayes Regression mache. Das geht mit dem Paket rstanarm recht einfach und man kann dann anhand der posterior-Verteilung ganz gut Vergleiche anstellen. Ich muss gerade arbeiten. Ich schicke Dir schonmal diesen Code und melde mich später nochmal mit mehr Erläuterungen.

CAVE: Das beruht auf einer MCMC-Simulation und wird daher etwas Rechenzeit brauchen. Lass das laufen und geh einen Kaffee kochen und auf die Toilette, nach zehn Minuten schreibt dieser Code Dir sechs Grafiken ins RStudio.

Code: Alles auswählen

d <- read.csv2("http://forum.r-statistik.de/download/file.php?id=1528")

library(rstanarm)
options(mc.cores = parallel::detectCores())

out1 <- list()
out2 <- list()
for(land in unique(d$land)){
  out1[[land]] <- stan_glm(x1~y1, data = d[d$land == land,])
  out2[[land]] <- stan_glm(x2~y2, data = d[d$land == land,])
}


par(mfrow=c(2,1))
for(i in unique(d$land)){
    hist(as.data.frame(out1[[i]])$y1, main = paste("y1 coefficent in land ", i), 
         xlab = expression(beta[1]))
    hist(as.data.frame(out2[[i]])$y2, main = paste("y2 coefficent in land ", i), 
         xlab = expression(beta[1]))
}
Bis später,
Bernhard
Rplot01.png
Rplot01.png (5.45 KiB) 580 mal betrachtet
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: Gleichheit von Koeffizienten ermitteln

Beitrag von bigben »

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:

Code: Alles auswählen

compare(4)
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
Rplot02.png (9.07 KiB) 568 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
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Antworten