Letzte Woche in Spielzimmer Andrew und ich haben an Code von gearbeitet der Oktober-Beitrag zum Thema „Äquivalentgewichte“Umfragegewichte, die einer Methode „äquivalent“ sind, die das Ergebnis Y modelliert (wie MRP).
Der Code vergleicht die Leistung von drei Arten von Umfragegewichten:
- Inverse-Response-Wahrscheinlichkeitsgewichte (IPW): W = 1/Ehat(R | X)
- kalibrierte Gewichte: W mit E(RWX) = E(X). Wenn X kategorisch ist, ist dies klassisch Poststratifizierung.
- äquivalente Gewichte: W mit E(RWY) = E(Ehat(Y | X, R=1)). Finden Sie Gewichtungen, sodass die Gewichtung der Stichprobe mit einer Schätzung übereinstimmt, die auf einem Modell für Ergebnis Y basiert (z. B. kleinste Quadrate oder mehrstufige Regression).
Mein ursprünglicher Code hat nur eine Simulation durchgeführt, hier machen wir 1000. Fortschritt!
library(lme4)
library(survey)
library(arm)
set.seed(1)
S = 1000
raw_means = rep(NA, S)
IPW_means = rep(NA, S)
cal_X_means = rep(NA, S)
equiv_LS_means = rep(NA, S)
equiv_MR_means = rep(NA, S)
true_means = rep(NA, S)
for (s in 1:S) {
Eintritt in die Simulationsschleife….
In jeder Simulation generieren wir Bevölkerungsdaten mit einem kategorischen X als Y | X ~ N(b_X, 6).
Okay = 100 N_per = 200 lev = issue(1:Okay) b_true = rnorm(Okay,0,3) names(b_true) = ranges(lev) pop = knowledge.body(X=rep(lev,every=N_per)) pop$Y = b_true(pop$X) + rnorm(nrow(pop),0,6) true_mean = imply(pop$Y) Npop = nrow(pop)
p_incl_k = exp(-0.8* (b_true-mean(b_true)))
p_incl_k = 0.02+0.18*p_incl_k/max(p_incl_k)
p_i = p_incl_k(pop$X) # true E(R|X) = f(b_X)
maintain=rbinom(Npop,1,p_i)==1
whereas (any(desk(pop$X(maintain)) ==0)) { maintain=rbinom(Npop,1,p_i)==1 }
samp=pop(maintain, )
Diese Stichprobe ist nicht vernachlässigbar, weil (Bayesianische Datenanalyse S.204):
Parameter des Datenerfassungsmechanismus unterscheiden sich nicht von den Parametern, die die Datenverteilung indizieren
Als nächstes erhalten wir Stichprobengrößen und geschätzte Einschlusswahrscheinlichkeiten für IPW:
nk=as.integer(desk(samp$X)) samp$nk=nk(samp$X) pop$nk=nk(pop$X) phat_incl_k=nk/N_per phat_i = phat_incl_k(pop$X) # Ehat(R|X) samp$w_hat=1/phat_i(maintain) pop$w_hat=1/phat_i
fit_ls=lm(Y~X,knowledge=samp) fit_mr=lmer(Y~1+ (1|X),knowledge=samp) samp$Yhat_LS=predict(fit_ls,newdata=samp) samp$Yhat_MR=predict(fit_mr,newdata=samp) pop$Yhat_LS=predict(fit_ls,newdata=pop) pop$Yhat_MR=predict(fit_mr,newdata=pop)
tot_X=as.vector(colSums(mannequin.matrix(~X-1,knowledge=pop))) names(tot_X) =colnames(mannequin.matrix(~X-1,knowledge=pop)) tot_Yhat_LS=sum(pop$Yhat_LS) tot_Yhat_MR=sum(pop$Yhat_MR) Npop=nrow(pop)
Dann verwenden wir Thomas Lumleys Umfragepaket:
d0=svydesign(ids=~1,weights=~1,knowledge=samp)
d_IPW=svydesign(ids=~1,weights=~w_hat,knowledge=samp)
pop_xtab=as.knowledge.body(desk(X=pop$X))
names(pop_xtab)(2) = "Freq"
d_X=postStratify(d0,~X,pop_xtab)
d_LS=calibrate(d0, ~1+Yhat_LS,
inhabitants=c(`(Intercept)`=Npop,Yhat_LS=tot_Yhat_LS),
calfun="linear")
d_MR=calibrate(d0, ~1+Yhat_MR,
inhabitants=c(`(Intercept)`=Npop,Yhat_MR=tot_Yhat_MR),
calfun="linear")
Und notieren Sie unsere Ergebnisse für diese eine Simulation:
raw_means(s) =coef(svymean(~Y,d0)) IPW_means(s) =coef(svymean(~Y,d_IPW)) cal_X_means(s) =coef(svymean(~Y,d_X)) equiv_LS_means(s) =coef(svymean(~Y,d_LS)) equiv_MR_means(s) =coef(svymean(~Y,d_MR)) true_means(s) =true_mean } # finish loop over simulations s
MSE_raw_mean = imply((raw_means - true_means)^2) MSE_IPW = imply((IPW_means - true_means)^2) MSE_cal_X = imply((cal_X_means - true_means)^2) MSE_equiv_LS = imply((equiv_LS_means - true_means)^2) MSE_equiv_MR = imply((equiv_MR_means - true_means)^2) MSEs = c( raw_mean=MSE_raw_mean, IPW=MSE_IPW, cal_X=MSE_cal_X, equiv_LS=MSE_equiv_LS, equiv_MR=MSE_equiv_MR ) spherical(MSEs, 3)
Die Ergebnisse zeigen, dass die Rohmittel mit Abstand am schlechtesten abschneiden. Wie erwartet sind IPW, Kalibrierung und äquivalente Gewichte der kleinsten Quadrate gleich (siehe Lumley 2010 Kapitel 7 und 9). Und die mehrstufige Regression funktioniert am besten (was ich nicht gesehen habe, als ich nur eine Simulation durchgeführt habe!). Würde MR besser abschneiden, wenn der Probenahmemechanismus vernachlässigbar wäre? Ich habe auch versucht, dem Modell Stichprobengrößen hinzuzufügen, wie Andrew vorgeschlagen hatte, aber die Ergebnisse wurden schlechter. Du bist dran!
raw_mean IPW cal_X equiv_LS equiv_MR 1.839 0.103 0.103 0.103. 0.099