Silakan pertimbangkan data ini:
dt.m <- structure(list(id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), occasion = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("g1", "g2"), class = "factor"), g = c(12, 8, 22, 10, 10, 6, 8, 4, 14, 6, 2, 22, 12, 7, 24, 14, 8, 4, 5, 6, 14, 5, 5, 16)), .Names = c("id", "occasion", "g"), row.names = c(NA, -24L), class = "data.frame")
Kami cocok dengan model komponen varians sederhana. Dalam R kita memiliki:
require(lme4)
fit.vc <- lmer( g ~ (1|id), data=dt.m )
Kemudian kami memproduksi plot ulat:
rr1 <- ranef(fit.vc, postVar = TRUE)
dotplot(rr1, scales = list(x = list(relation = 'free')))[["id"]]
Sekarang kami cocok dengan model yang sama di Stata. Pertama menulis ke format Stata dari R:
require(foreign)
write.dta(dt.m, "dt.m.dta")
Di Stata
use "dt.m.dta"
xtmixed g || id:, reml variance
Outputnya setuju dengan output R (tidak ditampilkan), dan kami berusaha untuk menghasilkan plot ulat yang sama:
predict u_plus_e, residuals
predict u, reffects
gen e = u_plus_e – u
predict u_se, reses
egen tag = tag(id)
sort u
gen u_rank = sum(tag)
serrbar u u_se u_rank if tag==1, scale(1.96) yline(0)
Clearty Stata menggunakan kesalahan standar yang berbeda untuk R. Bahkan Stata menggunakan 2,13 sedangkan R menggunakan 1,32.
Dari apa yang saya tahu, 1,32 dalam R berasal
> sqrt(attr(ranef(fit.vc, postVar = TRUE)[[1]], "postVar")[1, , ])
[1] 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977
meskipun saya tidak bisa mengatakan saya benar-benar mengerti apa yang dilakukan ini Bisakah seseorang menjelaskan?
Dan saya tidak tahu dari mana 2,13 dari Stata berasal, kecuali bahwa, jika saya mengubah metode estimasi ke kemungkinan maksimum:
xtmixed g || id:, ml variance
.... maka tampaknya menggunakan 1,32 sebagai kesalahan standar dan menghasilkan hasil yang sama dengan R ....
.... tetapi kemudian estimasi untuk varian efek acak tidak lagi setuju dengan R (35,04 vs 31,97).
Jadi sepertinya ada hubungannya dengan ML vs REML: Jika saya menjalankan REML di kedua sistem, output model setuju tetapi kesalahan standar yang digunakan dalam plot ulat tidak setuju, sedangkan jika saya menjalankan REML di R dan ML di Stata , plot ulat setuju, tetapi perkiraan model tidak.
Adakah yang bisa menjelaskan apa yang sedang terjadi?
sumber
[XT] xtmixed
dan / atau[XT] xtmixed postestimation
? Mereka merujuk pada Pinheiro dan Bates (2000), jadi setidaknya beberapa bagian matematika harus sama.Jawaban:
Menurut
[XT]
manual untuk Stata 11:Dari pertanyaan Anda, Anda telah mencoba REML di Stata dan R, dan ML di Stata dengan REML di R. Jika Anda mencoba ML di keduanya, Anda harus mendapatkan hasil yang sama di keduanya.
sumber