A followup the conversation that we were having. Basically, the idea that for a simulated, linear mixed effects model of the form \(Y = 0.5 + 0.5X_{1} + u_{0} + \epsilon_{ij}\) where \(u_{0} \sim N(0, 0.7)\) and \(\epsilon \sim N(0, 0.3)\) to get an intraclass correlation (ICC) of 0.7. ends up not showing anything even remotely heteroskedastic in the plot. Also, when you simulate data and test the residuals using the BreuschPagan test, the pvalue is not significant. HOWEVER the BreuschGodfrey of serial correlation *is* significant. Which low key took me by surprise because I assumed random intercept models are not necessarily autoregressive models of lag 1 (which is what the test in R shows... I think?)
Anyhoo, code:
Tests:
Anyhoo, code:
Code:
library(ggplot2)
library(lmtest)
set.seed(123)##para fig 4
library(lme4)
library(MASS)
Simulate < function(k, n, icc, beta, sigma=NA, sigmaXre=NA){
nk < n*k
nx < length(beta)1
Z < kronecker(diag(k), rep(1,n))
# X matrix
if (is.na(sigma[1])) sigma < diag(nx)
X < mvrnorm(nk, rep(0,nx), sigma)
# random effects of X
if (!is.na(sigmaXre[1])){
Xre < t(matrix(rnorm(k*nx,rep(1,nx),sigmaXre),nrow=nx))
Xre < cbind(rep(1,nk), X * (Z %*% Xre))
} else {
Xre < cbind(rep(1,nk), X)
}
X < cbind(rep(1,nk), X)
# create a factor to keep track of which "students" are in which "school"
group < as.factor(rep(1:k, each=n))
# generate taus and epsilons
tecov < diag(c(icc, rep(1icc,n)))
te < mvrnorm(k, rep(0,n+1), tecov)
epsilons < as.vector(t(te[,1]))
taus < te[,1]
# generate Y data
ran < Z %*% taus + epsilons
Y < Xre %*% beta + ran
output < list(Y, X[,2], group)
class(output) < "data.frame"
colnames(output) < c("Y", "X1", "group")
return(output)
}
# How many "schools"?
k < 100
# How many "students" per "school"?
n < 100
# Which ICC?
real.icc < 0.7
# Which fixed effects?
beta < c(0.5,0.5)
# Covariance matrices of the X variables (a positivedefinite symmetric matrix)
sigma < matrix(c(1), length(beta)1)
# Covariance matrix of the random effects of X (or vector of variances where zero means no random effect)
sigmaXre < matrix(c(rep(0,1)), length(beta)1) #interceptonly model
data1 < Simulate(k, n, real.icc, beta, sigma, sigmaXre)
Y <c(data1$Y)
X1 <c(data1$X1)
model < lm(Y~ X1)
Code:
> bptest(model)
studentized BreuschPagan test
data: model
BP = 0.14362, df = 1, pvalue = 0.7047
> ols_test_score(model)
Error in ols_test_score(model) : could not find function "ols_test_score"
> bgtest(model)
BreuschGodfrey test for serial correlation of order up to 1
data: model
LM test = 4729, df = 1, pvalue < 2.2e16
Attachments

331.5 KB Views: 4