[R] (no subject)

[Loai Mahmoud Awad Alzoubi]

Hi all

I am trying to generate a normal unbalanced data to estimate the coefficients of LM, LMM, GLM, and GLMM and their standard errors. Also, I am trying to estimate the variance components and their standard errors. Further, I am trying to use the likelihood ratio test to test H0: sigma^2_b = 0 (random effects variance component), and the t-test to test H0:mu=0 (intercept of the model Yij = mu + Bi + Eij).
I am using the following program
all.fix.coef.lme <- rep(NA,R) 
all.fix.coef.lm <- rep(NA,R) 
all.fix.coef.glm <- rep(NA,R) 
all.se.fix.coef <- rep(NA,R) 
all.fix.coef.glmm <- rep(NA,R) 
all.se.fix.coef.glmm <- rep(NA,R) 
all.varcomp.g <- rep(NA,R)
all.se.varcomp.g2 <- rep(NA,R)
all.se.varcomp.gD <- rep(NA,R)
all.se.fix.coef.adapt <- rep(NA,R)
vb<-rep(NA,R)
yb<-rep(NA,R)
va<-rep(NA,R)
mse<-rep(NA,R)
msa<-rep(NA,R)
maxH0<-rep(NA,R)
maxH1<-rep(NA,R)
tstat1 <- rep(NA,R)
tstat2 <- rep(NA,R)
tstatD<- rep(NA,R)
tstatad <- rep(NA,R)
tstatlme <- rep(NA,R)
tstatlm <- rep(NA,R)
sigmahatesq <- rep(NA,R)
sigmahatbsq <- rep(NA,R)
est.ICC <- rep(NA,R)
all.fix.coef.lm <- rep(NA,R) 
all.se.fix.coef.lm <- rep(NA,R)
all.se.fix.coef.glm <- rep(NA,R)
all.se.fix.coef.adapt.lm <- rep(NA,R) 
F<- rep(NA,R)
F1<- rep(NA,R)
lrt <- rep(NA,R)
lrtest <- rep(NA,R)
yb <- rep(NA,R) 
varyb <- rep(NA,R) 
var <- rep(NA,R) 
var1 <- rep(NA,R) 
D<- rep(NA,R)
Lrt <- function(k, R=250,m, c, stop=FALSE){
if(stop) browser()
set.seed(421)
options(digits = 3)
for(case in 1:3){       
if(case==1) c <- 1:3
if(case==2) c <- 5:15
if(case==3) c <- 25:150
for(r in(1:R)){ 
n <- c*m
       y <- rep(rnorm(m),each=c)+ rnorm(m*c)
       g <- as.factor(rep(1:m,each=c))
       test <- data.frame(cbind(g,y))
       #test$h <- 1:3
       test.lme1 <- lme (y~1, test, random = ~1|g)
       #test.lme2 <- lme (y~1, test)
       glm <- glm(y~1)
       test.lm <- lm (y~1)
       glmm <- glmmPQL (y~1, test, random = ~1|g, family=gaussian)
       all.se.fix.coef.lm[r] <- as.vector(sqrt(vcov(test.lm)))
       all.fix.coef.lme [r] <- as.vector(fixef(test.lme1))
       all.fix.coef.lm [r] <- as.vector(coef(test.lm))
       all.fix.coef.glm [r] <- as.vector(coef(glm))
       all.se.fix.coef [r] <- as.vector(sqrt(vcov(test.lme1)))
       all.fix.coef.glmm[r] <- as.vector(fixef(glmm)) 
       all.se.fix.coef.glmm[r] <- as.vector(sqrt(vcov(glmm)))
       lrtest[r]<- -2* (logLik(test.lme1, REML = TRUE))#- logLik(test.lme2, REML = TRUE))
if (lrtest[r]> 1.64) all.se.fix.coef.adapt[r]<-sqrt(vcov(test.lme1)) else all.se.fix.coef.adapt[r]<-sqrt(vcov(glm))
est.ICC[r] <- as.numeric(VarCorr(test.lme1)[1,1])/(        as.numeric(VarCorr(test.lme1)[1,1])+test.lme1$sigma^2) 
       all.varcomp.g [r] <-as.vector(as.numeric(VarCorr(test.lme1)[1,1]))
       if (test.lme1$apVar[1] == "Non-positive definite approximate variance-covariance") {
       all.se.varcomp.g2[r] <- NA
       tstat2[r] <- 0
 }
       else {
       all.se.varcomp.g2[r] <- sqrt(test.lme1$apVar[2,2])
       tstat2[r]<- all.varcomp.g[r]/all.se.varcomp.g2[r]
 }	
       all.se.varcomp.gD[r] <- sqrt(2/c*( c*as.numeric(VarCorr(test.lme1)[1,1]) + as.numeric(VarCorr(test.lme1)[2,1]))^2/(n+c)+(2/(c^2*(n-m+2))* as.numeric(VarCorr(test.lme1)[2,1])^2))
       tstatD[r] <- all.varcomp.g[r]/all.se.varcomp.gD[r]
       tstatlme[r]<- abs(all.fix.coef.lme[r]/all.se.fix.coef[r])
       tstatlm[r]<-abs(all.fix.coef.lm[r]/all.se.fix.coef.lm[r])  
       tstatad [r]<-abs(all.fix.coef.lme[r]/all.se.fix.coef.adapt[r])
 }
}
list(m, c,
            round(mean(all.se.fix.coef.adapt, na.rm = T), 3),
            round(mean(all.se.fix.coef.lm, na.rm = T ),3), 
            round(mean(all.se.fix.coef, na.rm = T ), 3),
            round(sqrt(var(all.fix.coef.lme, na.rm = T)),3), 
            round(mean(all.se.fix.coef.glmm), 3),
            #round(mean(all.fix.coef.glmm), 3),
            length(lrtest[lrtest>1.64])/2.50, 
            length(tstatlm[tstatlm>1.64])/2.5, 
            length(tstatlme[tstatlme>1.64] )/2.5, 
            length(tstatad[tstatad>1.64] )/2.5,
            round(mean(lrtest[lrtest>0]),3), 
            #round(var(lrtest[lrtest>0]), 3), 
            #length(lrtest[lrtest>0])/2.5, 
            length(tstatD[tstatD>1.64])/2.5, 
            length(lrtest[lrtest>1.64])/2.5,
            lrtest
)
  }
for (i in m) {
    for(j in c) {
     m = i
     c = j
A.i<- Lrt(0, ,m,c)
se.table <- rbind(se.table, A.i)
}
}
write.table(se.table, sep=”&”)

==========
But I couldn't find all values in the list, I just found 10 out of 14.

The results were as follows
"V1"&"V2"&"V3"&"V4"&"V5"&"V6"&"V7"&"V8"&"V9"&"V10"
2&2&0.489&0.457&0.514&0.532&10.8&34.4&28&32.8
2&5&0.323&0.302&0.303&0.355&6.8&41.2&33.6&39.2
2&10&0.235&0.219&0.224&0.256&6.4&56.4&47.2&51.6
2&50&0.108&0.099&0.099&0.118&6.8&99.6&94.8&97.6
2&100&0.074&0.071&0.073&0.082&3.2&100&99.2&99.2
5&2&0.316&0.303&0.307&0.334&12&40.8&33.2&38
5&5&0.201&0.197&0.203&0.215&4.4&69.2&62.4&68.4
5&10&0.147&0.14&0.142&0.156&7.2&86&80&82.8
5&50&0.066&0.063&0.067&0.07&6&100&100&100
5&100&0.047&0.045&0.044&0.05&6.4&100&100&100
10&2&0.224&0.219&0.228&0.232&8.4&62.4&59.6&62
10&5&0.144&0.14&0.139&0.15&8.4&88.4&84.4&86.8
10&10&0.103&0.099&0.1&0.107&8.8&98.8&97.6&98
10&50&0.046&0.045&0.046&0.048&7.2&100&100&100
10&100&0.033&0.032&0.033&0.034&8.8&100&100&100
50&2&0.101&0.1&0.1&0.102&8&98&98&98
50&5&0.064&0.063&0.066&0.065&8.4&100&100&100
50&10&0.045&0.045&0.046&0.046&8.4&100&100&100
50&50&0.02&0.02&0.02&0.021&7.2&100&100&100
50&100&0.014&0.014&0.014&0.015&9.2&100&100&100
100&2&0.071&0.071&0.074&0.072&6.4&100&100&100
100&5&0.045&0.045&0.045&0.046&11.6&100&100&100
100&10&0.032&0.032&0.029&0.032&8.4&100&100&100
100&50&0.014&0.014&0.014&0.014&8&100&100&100
100&100&0.01&0.01&0.01&0.01&6.8&100&100&100
2&2&0.492&0.46&0.523&0.538&10.8&33.2&28.4&31.2
2&5&0.323&0.303&0.308&0.36&6.4&44.8&37.6&42.4
2&10&0.237&0.22&0.228&0.261&7.2&64.8&56&61.2
2&50&0.119&0.099&0.123&0.131&13.2&98.4&91.2&93.2
2&100&0.09&0.071&0.106&0.101&17.2&100&97.2&97.2
5&2&0.319&0.305&0.31&0.336&13.2&44&38.4&41.6
5&5&0.203&0.198&0.206&0.217&5.6&73.6&68.4&72.4
5&10&0.15&0.141&0.144&0.159&11.2&90.8&82.8&86
5&50&0.072&0.063&0.077&0.078&20&100&100&100
5&100&0.058&0.045&0.064&0.062&34.8&100&100&100
10&2&0.226&0.221&0.23&0.235&9.2&68&65.6&67.6
10&5&0.145&0.141&0.144&0.153&8.4&92.8&89.6&92
10&10&0.105&0.1&0.104&0.11&12.4&99.6&98.8&98.8
10&50&0.052&0.045&0.056&0.055&28.8&100&100&100
10&100&0.042&0.032&0.047&0.044&54&100&100&100
50&2&0.101&0.1&0.103&0.103&11.2&98.8&98.4&98.8
50&5&0.065&0.064&0.067&0.067&14&100&100&100
50&10&0.047&0.045&0.048&0.048&21.6&100&100&100
50&50&0.024&0.02&0.024&0.024&71.6&100&100&100
50&100&0.02&0.014&0.02&0.02&98.4&100&100&100
100&2&0.071&0.071&0.075&0.072&8.4&100&100&100
100&5&0.046&0.045&0.047&0.047&17.2&100&100&100
100&10&0.033&0.032&0.032&0.033&24.4&100&100&100
100&50&0.017&0.014&0.017&0.017&90.4&100&100&100
100&100&0.014&0.01&0.015&0.014&100&100&100&100
2&2&0.507&0.467&0.547&0.551&12.8&38.8&32&35.2
2&5&0.331&0.306&0.337&0.374&6.8&48.8&41.6&46.4
2&10&0.257&0.222&0.26&0.285&12&70.8&56.4&64
2&50&0.167&0.1&0.196&0.178&34.8&94&77.2&79.2
2&100&0.143&0.072&0.181&0.151&42.4&97.2&84.8&86
5&2&0.327&0.31&0.321&0.345&15.6&50.4&42&48
5&5&0.21&0.2&0.223&0.228&9.2&77.2&71.6&72.8
5&10&0.163&0.143&0.166&0.177&21.2&92&80.8&85.6
5&50&0.108&0.064&0.116&0.112&63.2&100&99.6&99.6
5&100&0.102&0.046&0.112&0.103&83.6&100&100&100
10&2&0.234&0.227&0.238&0.244&12.4&72&67.6&70.8
10&5&0.154&0.144&0.16&0.163&18.8&93.2&90.4&92
10&10&0.118&0.102&0.12&0.124&32.8&100&98.4&98.4
10&50&0.08&0.046&0.086&0.082&87.6&100&99.6&99.6
10&100&0.076&0.032&0.079&0.076&97.6&100&100&100
50&2&0.104&0.102&0.109&0.107&17.6&99.2&99.2&99.2
50&5&0.07&0.065&0.072&0.072&44&100&100&100
50&10&0.054&0.046&0.057&0.055&75.6&100&100&100
50&50&0.038&0.021&0.037&0.038&100&100&100&100
50&100&0.036&0.015&0.036&0.036&100&100&100&100
100&2&0.073&0.072&0.077&0.074&19.6&100&100&100
100&5&0.05&0.046&0.051&0.05&56.8&100&100&100
100&10&0.039&0.032&0.038&0.039&94.8&100&100&100
100&50&0.027&0.015&0.027&0.027&100&100&100&100
100&100&0.025&0.01&0.026&0.025&100&100&100&100
2&2&0.677&0.573&0.877&0.758&23.2&49.6&41.6&45.2
2&5&0.574&0.362&0.757&0.627&36&61.2&43.6&48.4
2&10&0.596&0.266&0.698&0.619&55.2&72.8&45.2&47.6
2&50&0.581&0.12&0.752&0.586&78.4&86.4&47.6&48
2&100&0.545&0.084&0.707&0.549&82.8&92&47.6&48.4
5&2&0.497&0.415&0.507&0.523&48.4&61.6&44&49.2
5&5&0.432&0.259&0.472&0.441&77.2&74&54.8&55.6
5&10&0.441&0.188&0.455&0.443&92.8&84.4&57.6&57.6
5&50&0.417&0.083&0.441&0.417&100&94.4&58.8&58.8
5&100&0.414&0.059&0.444&0.414&99.6&96.8&58.4&58.4
10&2&0.374&0.311&0.379&0.385&68&76.8&67.6&68.8
10&5&0.337&0.194&0.351&0.338&97.6&91.6&74&74
10&10&0.324&0.137&0.306&0.324&99.2&96.8&80.4&80.4
10&50&0.306&0.061&0.315&0.306&100&99.6&80&80
10&100&0.307&0.043&0.308&0.307&100&99.6&78&78
50&2&0.173&0.141&0.18&0.173&100&99.6&99.6&99.6
50&5&0.154&0.089&0.151&0.154&100&100&100&100
50&10&0.149&0.063&0.155&0.149&100&100&100&100
50&50&0.143&0.028&0.138&0.143&100&100&100&100
50&100&0.143&0.02&0.146&0.143&100&100&100&100
100&2&0.122&0.099&0.121&0.122&100&100&100&100
100&5&0.109&0.063&0.104&0.109&100&100&100&100
100&10&0.105&0.045&0.105&0.105&100&100&100&100
100&50&0.101&0.02&0.099&0.101&100&100&100&100
100&100&0.101&0.014&0.106&0.101&100&100&100&100
100&c(25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, &0.102&0.028&0.102&0.132&0.101&100&70&20.4

So, would you please, help me doing that.

Loai ALzoubi
PhD candidate
UOW, Australia