2017-03-02 60 views
1

我正在使用多級模型嘗試描述縱向變化中的不同模式。當隨機效應完全相關時,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)) 
+0

總是可以貝葉斯。 –

+0

我正在考慮去貝葉斯,但只有有限的經驗與R2WinBUGS和R2jags。如果你能給我一個例子,它將不勝感激。 – Niek

+0

@MattTyers花了一些時間,但我嘗試了rjags包。儘管歡迎任何反饋 – Niek

回答

1

基於由@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