Я работаю с многоуровневыми моделями, пытаясь описать различные закономерности лонгитюдных изменений. Dingemanse et al (2010) описывают схему "разветвления", когда случайный эффекты полностью коррелированы. Однако я обнаружил, что аналогичная картина возникает, когда взаимосвязь между случайными эффектами нелинейна, а монотонно возрастает в наблюдаемом интервале. В этом случае случайные эффекты не полностью коррелированы, а описываются функцией. См. пример ниже для иллюстрации этого. Пример по-прежнему имеет высокую корреляцию точки пересечения и наклона (> 0,9), но можно получить корреляции ниже 0,7, сохраняя при этом идеальное соотношение точки пересечения и наклона.
Мой вопрос: есть ли способ включить идеальную (нелинейную) связь между случайными эффектами в многоуровневую модель, используя 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))
