在R中指定随机效果之间的非线性关系
问题描述:
我正在使用多级模型尝试描述纵向变化中的不同模式。当随机效应完全相关时,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))
答
基于由@MattTyers建议我使用rjags
贝叶斯方法一展身手。这是对随机效应之间的已知关系进行模拟数据的第一次尝试,但似乎产生了准确的估计值(好于nlme
模型)。我仍然对Gelman收敛诊断以及如何将此解决方案应用于实际数据感到担忧。但是,我想我会发布我的答案,以防有人正在解决同一问题。
# BAYESIAN ESTIMATE
library(ggplot2); library(reshape2)
# Set new dataset
set.seed(12345)
# New dataset to separate random and fixed
N<-100 # Number of respondents
int<-100 # Fixed effect intercept
U0<-rnorm(N,0,15) # Random effect intercept
slp_lin<-1 # Fixed effect linear slope
slp_qua<-.25 # Fixed effect quadratic slope
ID<- 1:100 # ID numbers
U1<-exp(U0/25)/7.5 # Random effect linear slope
U2<-exp(U0/25)/25 # Random effect quadratic slope
times<-15 # Max age
err <- matrix(rnorm(N*times,0,2),ncol = times, nrow = N) # Residual term
age <- 1:15 # Ages
# Create matrix of 'math' scores using model
math<-matrix(NA,ncol = times, nrow = N)
for(i in ID){
for(j in age){
math[i,j] <- (int + U0[i]) +
(slp_lin + U1[i])*age[j] +
(slp_qua + U2[i])*(age[j]^2) +
err[i,j]
}}
# Melt dataframe and plot scores
e.long<-melt(math)
names(e.long) <- c("ID","age","math")
ggplot(e.long,aes(x= age, y= math)) + geom_line(aes(group = ID))
# Create dataframe for rjags
dat<-list(math=as.numeric(e.long$math),
age=as.numeric(e.long$age),
childnum=e.long$ID,
n=length(e.long$math),
nkids=length(unique(e.long$ID)))
lapply(dat , summary)
library(rjags)
# Model with uninformative priors
model_rnk<-"
model{
#Model, fixed effect age and random intercept-slope connected
for(i in 1:n)
{
math[i]~dnorm(mu[i], sigm.inv)
mu[i]<-(b[1] + u[childnum[i],1]) + (b[2]+ u[childnum[i],2]) * age[i] +
(b[3]+ u[childnum[i],3]) * (age[i]^2)
}
#Random slopes
for (j in 1:nkids)
{
u[j, 1] ~ dnorm(0, tau.a)
u[j, 2] <- exp(u[j,1]/25)/7.5
u[j, 3] <- exp(u[j,1]/25)/25
}
#Priors on fixed intercept and slope priors
b[1] ~ dnorm(0.0, 1.0E-5)
b[2] ~ dnorm(0.0, 1.0E-5)
b[3] ~ dnorm(0.0, 1.0E-5)
# Residual variance priors
sigm.inv ~ dgamma(1.5, 0.001)# precision of math[i]
sigm<- pow(sigm.inv, -1/2) # standard deviation
# Varying intercepts, varying slopes priors
tau.a ~ dgamma(1.5, 0.001)
sigma.a<-pow(tau.a, -1/2)
}"
#Initialize the model
mod_rnk<-jags.model(file=textConnection(model_rnk), data=dat, n.chains=2)
#burn in
update(mod_rnk, 5000)
#collect samples of the parameters
samps_rnk<-coda.samples(mod_rnk, variable.names=c("b","sigma.a", "sigm"), n.iter=5000, n.thin=1)
#Numerical summary of each parameter:
summary(samps_rnk)
gelman.diag(samps_rnk, multivariate = F)
# nlme model
library(nlme)
Stab_rnk2<-lme(math ~ age + I(age^2),
random= list(ID = ~age + I(age^2)),
data=e.long,method="ML",na.action=na.exclude)
summary(Stab_rnk2)
结果看起来非常接近人口值
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
b[1] 100.7409 0.575414 5.754e-03 0.1065523
b[2] 0.9843 0.052248 5.225e-04 0.0064052
b[3] 0.2512 0.003144 3.144e-05 0.0003500
sigm 1.9963 0.037548 3.755e-04 0.0004056
sigma.a 16.9322 1.183173 1.183e-02 0.0121340
而且NLME估计是差远了(与随机拦截的除外)
Random effects:
Formula: ~age + I(age^2) | ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 16.73626521 (Intr) age
age 0.13152688 0.890
I(age^2) 0.03752701 0.924 0.918
Residual 1.99346015
Fixed effects: math ~ age + I(age^2)
Value Std.Error DF t-value p-value
(Intercept) 103.85896 1.6847051 1398 61.64816 0
age 1.15741 0.0527874 1398 21.92586 0
I(age^2) 0.30809 0.0048747 1398 63.20204 0
总是可以贝叶斯。 –
我正在考虑去贝叶斯,但只有有限的经验与R2WinBUGS和R2jags。如果你能给我一个例子,它将不胜感激。 – Niek
@MattTyers花了一些时间,但我尝试了rjags包。尽管欢迎任何反馈 – Niek