Mengapa saya mendapatkan hasil yang sangat berbeda untuk poli (mentah = T) vs poli ()?

10

Saya ingin memodelkan dua variabel waktu yang berbeda, beberapa di antaranya sangat collinear dalam data saya (usia + kohort = periode). Melakukan ini saya mengalami beberapa masalah dengan lmerdan dan interaksi poly(), tetapi mungkin tidak terbatas lmer, saya mendapat hasil yang sama dengan nlmeIIRC.

Jelas, pemahaman saya tentang apa fungsi poly () tidak ada. Saya mengerti apa yang poly(x,d,raw=T)dilakukan dan saya pikir tanpanya raw=Tmembuat polinomial ortogonal (saya tidak bisa mengatakan saya benar-benar mengerti apa artinya), yang membuat pemasangan lebih mudah, tetapi tidak membiarkan Anda menafsirkan koefisien secara langsung.
Saya membaca bahwa karena saya menggunakan fungsi prediksi, prediksi harus sama.

Tetapi mereka tidak, bahkan ketika model bertemu secara normal. Saya menggunakan variabel terpusat dan saya pertama kali berpikir bahwa mungkin polinomial ortogonal mengarah ke korelasi efek tetap yang lebih tinggi dengan istilah interaksi collinear, tetapi tampaknya sebanding. Saya telah menyisipkan dua ringkasan model di sini .

Plot-plot ini diharapkan menggambarkan perbedaan yang ada. Saya menggunakan fungsi prediksi yang hanya tersedia di dev. versi lme4 (mendengarnya di sini ), tetapi efek tetapnya sama dalam versi CRAN (dan mereka juga kelihatan tidak nyaman, misalnya ~ 5 untuk interaksi ketika DV saya memiliki kisaran 0-4).

Panggilan singkatnya adalah

cohort2_age =lmer(churchattendance ~ 
poly(cohort_c,2,raw=T) * age_c + 
ctd_c + dropoutalive + obs_c + (1+ age_c |PERSNR), data=long.kg)

Prediksi ini adalah efek tetap saja, pada data palsu (semua prediktor lain = 0) di mana saya menandai rentang yang ada dalam data asli sebagai ekstrapolasi = F.

predict(cohort2_age,REform=NA,newdata=cohort.moderates.age)

Saya dapat memberikan lebih banyak konteks jika perlu (saya tidak berhasil menghasilkan contoh yang dapat direproduksi dengan mudah, tetapi tentu saja dapat berusaha lebih keras), tetapi saya pikir ini adalah permintaan yang lebih mendasar: jelaskan poly()fungsinya kepada saya, tolong.

Polinomial mentah

Polinomial mentah

Polinomial ortogonal (terpotong, tidak terpotong di Imgur )

Polinomial ortogonal

Ruben
sumber

Jawaban:

10

Saya pikir ini adalah bug dalam fungsi prediksi (dan karenanya kesalahan saya), yang notabene nlme tidak berbagi. ( Sunting : harus diperbaiki di versi R-forge terbaru lme4.) Lihat di bawah untuk contoh ...

Saya pikir pemahaman Anda tentang polinomial ortogonal mungkin baik-baik saja. Hal rumit yang perlu Anda ketahui tentang mereka jika Anda mencoba menulis metode prediksi untuk kelas model adalah bahwa dasar untuk polinomial ortogonal ditentukan berdasarkan set data yang diberikan, jadi jika Anda naif (seperti saya lakukan! ) gunakan model.matrixuntuk mencoba menghasilkan matriks desain untuk set data baru, Anda mendapatkan basis baru - yang tidak lagi masuk akal dengan parameter lama. Sampai saya memperbaiki ini, saya mungkin perlu menjebak yang memberitahu orang-orang yang predicttidak bekerja dengan basis polinomial ortogonal (atau basis spline, yang memiliki properti yang sama).

d <- expand.grid(x=seq(0,1,length=50),f=LETTERS[1:10])
set.seed(1001)
u.int <- rnorm(10,sd=0.5)
u.slope <- rnorm(10,sd=0.2)
u.quad <- rnorm(10,sd=0.1)
d <- transform(d,
               ypred = (1+u.int[f])+
               (2+u.slope[f])*x-
               (1+u.quad[f])*x^2)
d$y <- rnorm(nrow(d),mean=d$ypred,sd=0.2)
ggplot(d,aes(x=x,y=y,colour=f))+geom_line()+
    geom_line(aes(y=ypred),linetype=2)

library(lme4)
fm1 <- lmer(y~poly(x,2,raw=TRUE)+(1|f)+(0+x|f)+(0+I(x^2)|f),
            data=d)


fm2 <- lmer(y~poly(x,2)+(1|f)+(0+x|f)+(0+I(x^2)|f),
            data=d)
newdat <- data.frame(x=unique(d$x))
plot(predict(fm1,newdata=newdat,REform=NA))
lines(predict(fm2,newdata=newdat,REform=NA),col=2)
detach("package:lme4")

library(nlme)
fm3 <- lme(y~poly(x,2,raw=TRUE),
           random=list(~1|f,~0+x|f,~0+I(x^2)|f),
            data=d)
VarCorr(fm3)

fm4 <- lme(y~poly(x,2),
           random=list(~1|f,~0+x|f,~0+I(x^2)|f),
            data=d)

newdat <- data.frame(x=unique(d$x))
lines(predict(fm3,newdata=newdat,level=0),col=4)
lines(predict(fm4,newdata=newdat,level=0),col=5)
Ben Bolker
sumber
Terima kasih, itu meyakinkan. Untuk mengulangi: Saya membaca bahwa Anda tidak dapat mengambil efek tetap polinomial ortogonal pada nilai nominal, tetapi kadang-kadang mereka tampak sangat besar. Misalnya jika saya menjalankan interaksi dua polinomial kubik, saya mendapatkan efek tetap untuk polinomial dan interaksinya dalam kisaran dari -22 hingga -127400. Itu sepertinya jauh bagi saya, terutama mengingat bahwa semua efek tetap adalah negatif. Akankah fungsi prediksi yang direvisi masuk akal dari efek-efek tetap ini atau apakah model-model tersebut secara palsu bersatu atau apakah ada yang salah pada akhirnya?
Ruben
Sekali lagi, saya curiga (tapi jelas tidak tahu pasti) bahwa semuanya baik-baik saja. Orth. polinomial baik untuk stabilitas numerik dan pengujian hipotesis, tetapi (saat Anda mengetahuinya) nilai parameter aktual bisa lebih sulit untuk ditafsirkan. Versi lme4-devel saat ini (saya baru saja memposting versi yang harus lulus tes, mungkin butuh ~ 24 jam untuk membangun kembali di r-forge, kecuali Anda dapat membuat sendiri dari SVN) harus memberi Anda prediksi yang cocok antara polinomial mentah / orto. Alternatifnya adalah memusatkan dan mengukur prediktor berkesinambungan à la Schielzeth 2010 Metode dalam Ekologi & Evolusi ...
Ben Bolker
Ya, kedua polinomial sangat setuju sekarang. Terima kasih banyak! Saya telah meningkatkan dan memusatkan prediktor saya, tetapi beberapa model tidak cocok dengan polinomial mentah.
Ruben