Dear All,
I would be very grateful for your help concerning the following question:
Below mentioned programme is available on net to generate longitudinal
data. Usually we get almost same parameter estimates as used to generate
the data. The problem here is I am not able to get it for data used here,
despite increasing sample size and number of simulations. Is it normal to
expect this type of behavior from mixed effect models? The problem is this
code running fine for data(#data) provided by the author. I am not good
with programming, can anyone comments on the code or
possible explanation for this behavior?
Any help would be greatly appreciated.
sim.lmer <- function(n,p,error,B0,varb0,B1,varb1,cor)
{
#start function
###########################################
########### INPUTS ####################
###########################################
#n is number of subjects
#p is number of time points
#error is Residuals of measurement
#B0 is fixed intercept effect (average group intercept)
#B1 is fixed slope effect (average group slope)
#varb0 is the variance of individual intercepts
#varb1 is the variance of individual slopes
#cor is the correlation between ind. intercepts and ind. slopes
################# ASSUMPTIONS#######################
# Random effects have a bivariate normal distribution with zero mean
# Random errors have a normal distribution with zero mean
# The responses are generated based on a linear mixed effect regression
model
# Y = B0 + B1*Weeks + b0 + b1 + e
#####################################################
#####################################################
#####################################################
# Start function
cov <- cor*sqrt(varb0*varb1) #correlation between random effects
d <- matrix(c(varb0,cov,cov,varb1),nrow=2,ncol=2) #var-cov matrix of random
effects
ind <- mvrnorm(n,c(0,0),d) #generate bivariate normal random
effects
b0 <- ind[,1] # individual intercepts' deviation from fixed intercept
b1 <- ind[,2] # individual slopes' deviation from fixed slope
d2 <- (error^2)*diag(p) # var-cov matrix of error terms at time points
err <- mvrnorm(n,rep(0,p),d2) # generate multivariate normal error terms
with zero mean
ind.slo <- B1+b1
ind.int <- B0+b0
rand.eff <- cbind(ind.slo,ind.int)
##########################################
data <- matrix(nrow=n,ncol=p)
for(i in 1:n) {
for(k in 1:p) {
data[i,k]= B0 + B1*(k-1) + (b0[i] + b1[i]*(k-1)) }}
data2 <- data + err
##############################################
data2 <- as.data.frame(data2)
names <- c()
for(i in 1:p) names[i]=paste("Score",i,sep="")
colnames(data2) <- names ### Add column names to data
mynames <- colnames(data2)
data2$ID <- 1:n
d <- reshape(data2,varying=mynames,idvar="ID",
v.names="Score",timevar="Weeks",times=1:p,direction="long")
list(data=d,rand.eff=rand.eff)
} ### End of Function
require(MASS)
require(plyr)
require(lme4)
data <-
sim.lmer(n=1000,p=6,error=10,B0=38.93,varb0=30.92,B1=-2.31,varb1=0.56,cor=-0.7)
#data <-
sim.lmer(n=1000,p=10,error=5,B0=40,varb0=1050,B1=.2,varb1=.4,cor=.7)
d <- data$data
LMER.1 <- lmer(Score ~ Weeks + (1 + Weeks | ID), data = d, REML=F)
summary(LMER.1)
--
Thanks & Regards,
Kamal Kishore
PhD Student
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
