2015-11-15 17 views
5

Juliaのlmer関数と同じ結果を、lme4ライブラリから取得したいと考えています。 Rを使用してmtcarsデータセットを作成してくださいR(lme4)からJuliaまでの混合エフェクトのモデル式を書き直す

library(lme4) 
data<-mtcars  
data$vs<-as.factor(data$vs) 
data$am<-as.factor(data$am) 
data$gear<-as.factor(data$gear) 
str(data) 
model <- lmer(mpg ~ cyl:gear + hp:am + (1|gear:am), data = data) 

私は、しかし、私はlmm()を使用してlmer()関数の最初の引数から式を書き換える方法がわからない、同じ結果を実行することが可能であるべきでジュリアのためのMixedModelsパッケージからlmm()機能を発見しました。特にinteraction(:)演算子。

私は短い例で答えに感謝します。

+0

これはhttps://github.com/dmbates/MixedModels.jlを助けるかもしれない – MLavoie

+2

(あなたが唯一の対話をしたい場合)の相互作用演算子は、ちょうど '*' '(これはあなたの主効果との相互作用を示します)や'&されます。 – mmagnuski

答えて

2

ここのRとジュリアモデルの対応は正確ではないようです。さまざまな数値アルゴリズムの問​​題もあります。しかし、次のようにMixedModelsを使用して同じモデルを再作成する私の試みは、次のとおりです。

using RDatasets 
using MixedModels 

mtcars = dataset("datasets","mtcars") 
mtcars[:AM] = PooledDataArray(mtcars[:AM]) 
mtcars[:Gear] = PooledDataArray(mtcars[:Gear]) 
mtcars[:GearAM] = PooledDataArray(collect(zip(mtcars[:Gear],mtcars[:AM]))) 
m = fit!(lmm(MPG ~ 1 + Gear + AM + Cyl + HP + (1|GearAM),mtcars)) 

が手動で混合効果列の作成klunkyある - おそらくより良い方法があります。 RとJuliaの間の係数の命名の違いに注意してください。両方に6つの固定効果係数があります。

ジュリアのソリューションは私のマシン上では違っているようですが、より良い対数尤度が得られます。ランダム効果は、その変数が既に固定効果(ギアとAMの間の依存性のみを考慮している)に存在するため、弱いと予想され、わずか32のデータ点しか存在しない。

これがうまくいくと思いますが、理解を深めることができれば、別の回答やコメントに追加するとよいでしょう。

関連する問題