AIC, anova error: Model tidak semua cocok dengan jumlah pengamatan yang sama, model tidak semua dipasang dengan ukuran dataset yang sama

Saya punya model seperti ini:

require(nlme)

set.seed(123)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, each = k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)

m1 <- lm(y ~ x)
summary(m1)

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

Sekarang saya mencoba menilai apakah efek acak harus ada dalam model. Jadi saya membandingkan model menggunakan AIC atau anova, dan saya mendapatkan kesalahan berikut:

> AIC(m1, m2, m3)
   df       AIC
m1  3 1771.4696
m2  7 -209.1825
m3  4 -154.0245
Warning message:
In AIC.default(m1, m2, m3) :
  models are not all fitted to the same number of observations  
> anova(m2, m3)
Error in anova.lmlist(object, ...) : 
  models were not all fitted to the same size of dataset

Seperti yang Anda lihat, dalam kedua kasus saya menggunakan dataset yang sama. Saya telah menemukan dua solusi, tetapi saya tidak menganggapnya memuaskan:

1. Menambah method = "ML"panggilan lme () - tidak yakin apakah ada baiknya mengubah metode.
2. Menggunakan lmer()sebagai gantinya. Anehnya, ini bekerja, meskipun faktanya lmer () menggunakan metode REML. Namun saya tidak suka solusi ini karena lmer()tidak menunjukkan nilai p untuk koefisien - saya lebih suka menggunakan yang lebih tua lme().

Apakah Anda tahu apakah ini bug atau tidak dan bagaimana kita bisa mengatasinya?

Pencarian cepat menunjukkan bahwa itu mungkin (walaupun saya harus mengakui bahwa saya pikir itu bukan) dan itu bukan bug ... hanya kasus lain di mana metode dalam R disembunyikan dan menghasilkan hal-hal yang tampaknya 'tidak terduga' ', tetapi kerumunan RTFM mengatakan, "Itu ada di dokumentasi." Bagaimanapun ... solusi Anda adalah anovadengan menggunakan lmeargumen pertama dan lmmodel sebagai argumen kedua (dan ketiga jika Anda suka). Jika ini tampak aneh, itu karena itu agak aneh. Alasannya adalah bahwa ketika Anda menelepon anova, anova.lmemetode ini dipanggil hanya jika argumen pertama adalah lmeobjek. Kalau tidak, itu panggilan anova.lm(yang pada gilirannya panggilan anova.lmlist; jika Anda menggali anova.lm, Anda akan melihat alasannya). Untuk perincian tentang bagaimana Anda ingin meneleponanovadalam hal ini, tarik bantuan untuk anova.lme. Anda akan melihat bahwa Anda dapat membandingkan model lain dengan lmemodel, tetapi mereka harus berada di posisi selain argumen pertama. Rupanya itu juga mungkin untuk digunakan anovapada model yang cocok menggunakan glsfungsi tanpa terlalu peduli tentang urutan argumen model. Tapi saya tidak tahu cukup detail untuk menentukan apakah itu ide yang baik atau tidak, atau apa sebenarnya implikasinya (tampaknya mungkin baik-baik saja, tetapi panggilan Anda). Dari tautan yang dibandingkan lmdengan lmetampaknya didokumentasikan dengan baik dan dikutip sebagai metode, jadi saya akan salah arah itu, kan saya.

Semoga berhasil.

Ini aneh tentunya. Sebagai pemikiran pertama: ketika melakukan perbandingan model di mana model memiliki struktur efek tetap yang berbeda ( m2dan m3misalnya), yang terbaik bagi kami adalah karena akan "berubah" . (Ini akan melipatgandakannya dengan , di mana ) Sangat menarik bahwa ia bekerja menggunakan yang membuat saya percaya itu mungkin bukan bug. Sepertinya hampir seperti menegakkan "praktik yang baik".y k k X = $ML$REML$y$$k$$kX= 0$method="ML"

Karena itu, mari kita lihat di bawah tenda:

 methods(AIC)  
 getAnywhere('AIC.default')

 A single object matching ‘AIC.default’ was found
 It was found in the following places
   registered S3 method for AIC from namespace stats
   namespace:stats with value

 function (object, ..., k = 2) 
 {
     ll <- if ("stats4" %in% loadedNamespaces()) 
         stats4:::logLik
     else logLik
     if (!missing(...)) {
         lls <- lapply(list(object, ...), ll)
         vals <- sapply(lls, function(el) {
             no <- attr(el, "nobs") #THIS IS THE ISSUE!
             c(as.numeric(el), attr(el, "df"), if (is.null(no)) NA_integer_ else no)
         })
         val <- data.frame(df = vals[2L, ], ll = vals[1L, ])
         nos <- na.omit(vals[3L, ])
         if (length(nos) && any(nos != nos[1L])) 
             warning("models are not all fitted to the same number of observations")
         val <- data.frame(df = val$df, AIC = -2 * val$ll + k * val$df)
             Call <- match.call()
             Call$k <- NULL
         row.names(val) <- as.character(Call[-1L])
         val
     }
     else {
         lls <- ll(object)
         -2 * as.numeric(lls) + k * attr(lls, "df")
     }     
 }

di mana dalam kasus Anda, Anda dapat melihat bahwa:

  lls <- lapply(list(m2,m3), stats4::logLik)
  attr(lls[[1]], "nobs")
  #[1] 500
  attr(lls[[2]], "nobs")
  #[1] 498

dan jelas logLikmelakukan sesuatu (mungkin?) tak terduga ...? tidak ada, tidak benar-benar, jika Anda melihat dokumentasi logLik, ?logLik, Anda akan melihatnya secara eksplisit dinyatakan:

 There may be other attributes depending on the method used: see
 the appropriate documentation.  One that is used by several
 methods is ‘"nobs"’, the number of observations used in estimation
 (after the restrictions if ‘REML = TRUE’)

yang membawa kami kembali ke titik awal kami, Anda harus menggunakan ML.

Untuk menggunakan ungkapan umum dalam CS: "Ini bukan bug; ini fitur (nyata)!"

EDIT : (Hanya untuk menjawab komentar Anda :) Asumsikan Anda cocok dengan model lain menggunakan lmerwaktu ini:

m3lmer <- lmer(y ~ x + 1|cat)

dan Anda melakukan hal berikut:

lls <- lapply(list(m2,m3, m3lmer), stats4::logLik)
attr(lls[[3]], "nobs")
#[1] 500
 attr(lls[[2]], "nobs")
#[1] 498

Yang tampak seperti perbedaan yang jelas antara keduanya tetapi sebenarnya tidak seperti yang dijelaskan Gavin. Hanya untuk menyatakan yang sudah jelas:

 attr( logLik(lme(y ~ x, random = ~ 1|cat, na.action = na.omit, method="ML")),
 "nobs")
#[1] 500

Ada alasan bagus mengapa ini terjadi dalam hal metodologi saya pikir. lmetidak mencoba untuk memahami regresi LME untuk Anda lmersaat melakukan perbandingan model, itu kembali segera ke hasil ML itu. Saya pikir tidak ada bug dalam hal ini lmedan lmerhanya alasan yang berbeda di balik setiap paket.

Lihat juga komentar Gavin Simposon tentang penjelasan yang lebih mendalam tentang apa yang terjadi anova()(Hal yang sama terjadi secara praktis dengan AIC)