我正在使用多級模型嘗試描述縱向變化中的不同模式。當隨機效應完全相關時,Dingemanse et al (2010)描述了「扇出」模式。然而,我發現當隨機效應之間的關係是非線性的但在觀察到的時間間隔內單調遞增時會出現類似的模式。在這種情況下,隨機效應並不完全相關,而是由函數描述。 請參閱下面的示例以獲取此示例。這個例子仍然具有很高的截距 - 斜率相關性(> .9),但是可以得到低於.7的相關性,同時仍然保持完美的截距 - 斜率關係。在R中指定隨機效果之間的非線性關係
我的問題:是否有辦法在多級模型中使用nlme
或其他一些R包包含隨機效應之間的完美(非線性)關係? MLwiN有一種方法來約束斜率截距協方差,這將是一個開始,但不足以包含非線性關係。到目前爲止,我一直無法找到nlme
的解決方案,但也許你知道一些其他可以在模型中包含此功能的軟件包。
道歉的草率編碼。我希望我的問題足夠清楚,但是如果需要澄清,請告訴我。任何幫助或替代解決方案,不勝感激。
set.seed(123456)
# Change function, quadratic
# Yit = B0ij + B1ij*time + B2ij*time^2
chn <- function(int, slp, slp2, time){
score<-int + slp * time+ slp2 * time^2
return(score)
}
# Set N, random intercept, time and ID
N<-100
start<-rnorm(N,100,15) # Random intercept
time<- matrix(1:15,ncol = 15, nrow = 100,byrow = T) # Time, balanced panel data
ID<-1:N # ID variable
# Random intercept, linear slope: exp(intercept/25)/75, quadratic slope: exp(intercept/25)/250
score3<-matrix(NA,ncol = ncol(time), nrow = N)
for(x in ID){
score3[x,]<-chn(start[x],exp(start[x]/25)/75,exp(start[x]/25)/250,time[x,])
}
#Create dataframe
df<- data.frame(ID,score3)
df2<- melt(df,id = 'ID')
df2$variable<-as.vector(time)
# plot
ggplot(df2, aes(x= variable, y= value)) + geom_line(aes(group = ID)) +
geom_smooth(method = "lm", formula = y ~ x + I(x^2), se =F, size = 2, col ="red")
# Add noise and estimate model
df2$value2<-df2$value + rnorm(N*ncol(time),0,2)
# Random intercept
mod1<-lme(value2 ~ variable + I(variable^2),
random= list(ID = ~1),
data=df2,method="ML",na.action=na.exclude)
summary(mod1)
# Random slopes
mod2<-update(mod1,.~.,random= list(ID = ~variable + I(variable^2)))
summary(mod2)
pairs(ranef(mod2))
總是可以貝葉斯。 –
我正在考慮去貝葉斯,但只有有限的經驗與R2WinBUGS和R2jags。如果你能給我一個例子,它將不勝感激。 – Niek
@MattTyers花了一些時間,但我嘗試了rjags包。儘管歡迎任何反饋 – Niek