Dear all, <r-sig-mixed-models@r-project.org>

I have a question about my simulation study of multilevel data. Concretely,
I want to simulate a model with one micro and one macro variable
y=Gamma00-Gamma01xij+Gamma10zj+Gamma11xijzj+u0j+u1j+rij. This are my
'estimates':

Gamma00=1.0 #average intercept across the level-2 units

Gamma10=0.3 #average regression slope across the level-2 units

Gamma01=0.3

Gamma11=0.3

 VAR.mikro=0.5

SD.mikro=sqrt(0.5)

VAR.interc1=0.05555556

SD.interc1=0.2357023

VAR.slope1=0.05555556

SD.slope1=0.2357023

cov(u0j,u1j) <- 0

 x <- rnorm (n, 0, 1)

z <- rnorm (J, 0, 1)

I've programmed the code in R, however I get wrong results, that is
correlation between intercept and slope is +1 or -1 and when I want to
run simulations, I get the message that the sigma is not positive
definite.  Below is the R code:



library(MASS)

library(lme4)

library(arm)



set.seed(1)


*#fixed effect*

npred<-2 #the number of fixed predictors (not including the intercept)

randsize<-2 #the number of unique elements in the variance matrix at level 2

beta<-c(1,0.3,0.3,0.3)

betasize<-length(beta) #length

effectbeta<-abs(beta) #effect size

sgnbeta<-sign(beta) #sign



*#random slope model*

randcolumn<-c(2.000000)



# store the mean and variances of the predictors (at the first and second
levels)

meanpred<-c(0,0.000000,0.000000)

varpred<-c(0,1.000000,0.000000)

varpred2<-c(0,1.000000,1.000000)

#define the variances of the residuals at the second and first levels,
respectively.

sigma2u<-matrix(c(0.05,0.0,0.0,0.05),randsize,randsize)

sigmae<-sqrt(0.500000)



##-----------------        Inputs for model fitting
-----------------##



fixname<-c("x0","x","z","xz") #The predictors are specified by the variable
fixname

fixform<-"1+x+z+ x:z"

randform<-"(1+x|group)"

expression<-paste(c(fixform,randform),collapse="+") #combining the form of
the fixed and random parts (fixform, randform, expression)

modelformula<-formula(paste("y ~",expression))

data<-vector("list",2+length(fixname))

names(data)<-c("group","y",fixname)

* *

*n1<-5*

*n2<-30*

length=n1*n2

x<-matrix(1,length,betasize)

z<-matrix(1,length,randsize)

group<-rep(c(1:n2),each=n1)



#######------------ To set up X matrix --------------########

* *

#mikro spremenljivka na mikro in makro ravni

micpred<-rnorm(length,meanpred[2],sqrt(varpred[2]))

macpred<-rnorm(n2,0,sqrt(varpred2[2]))

macpred<-rep(macpred,each=n1)

x[,2]<-micpred+macpred



#makro spremenljivka na mikro=0! in makro ravni, namen za generalizacijo!

micpred<-rnorm(length,meanpred[3],sqrt(varpred[3]))

macpred<-rnorm(n2,0,sqrt(varpred2[3]))

macpred<-rep(macpred,each=n1)

x[,3]<-micpred+macpred



#cross-level interaction

x[,4]<- x[,2]* x[,3]

#######------------       To set up Z matrix          --------------########



for(j in 2:dim(z)[2])

z[,j]<-x[,randcolumn[j-1]]



#####-----------------------------------------------------------------------#####




e<-rnorm(length,0,sigmae)

u<-mvrnorm(n2,rep(0,randsize),sigma2u)

fixpart<-x%*%beta

randpart<-rowSums(z*u[group,])

y<-fixpart+randpart+e



##------------------- Inputs for model fitting ---------------##



data$group<-as.factor(group)

data$y<-y

data$x0<-x[,1]

data$x<-x[,2]

data$z<-x[,3]

data$xz<-x[,4]



###~~~~~~~~~~ Fitting the model using lmer funtion ~~~~~~~~~~###



(fitmodel <- lmer(modelformula,data))


sim(fitmodel)

If anyone knows where the problem lies I will be very happy. Thank you.

	[[alternative HTML version deleted]]


