Dear all,
My code is presented as the following.
library(MASS)
library(rmutil)
library(repeated)
library(lme4)
library(arm)
#install.packages("Zelig", repos = "http://gking.harvard.edu";)
library(Zelig)
rm(list = ls())
beta0<-2.5
beta1<--0.3
sigs2<-0.5
I<-4
#numberpatients<-c(40,100,160,200,400,600)
numberpatients<-c(1000)
#numberpatients<-3
times<-1
MSE_eBLP_quad<-NULL
for (i3 in 1:length(numberpatients)) {
n<-numberpatients[i3]
b0<-NULL
b1<-NULL
sigma2<-NULL
eBLP_glimmix_quad<-NULL
for (kk in 1:times) {
s<-rnorm(I,0,sqrt(sigs2))
y<-NULL
u<-NULL
center<-NULL
for (i in 1:I) {
for (k in 1:n) {
u<-c(u,0)
center<-c(center,i)
lam<-exp(beta0+s[i])
obs<-rpois(1,lam)
y<-c(y,obs)
}
for (k in 1:n) {
u<-c(u,1)
center<-c(center,i)
lam<-exp(beta0+beta1+s[i])
obs<-rpois(1,lam)
y<-c(y,obs)
}
}
df<-data.frame(y,u,center)
##method1##############
#dan<-glmm(y~u,family=poisson,nest=center,points=3,data=df,print.level=2)
#names(dan)
#b0<-dan$coef[1]
#b1<-dan$coef[2]
#sigma2<-dan$coef[3]
yi1<-NULL
yi2<-NULL
for (i2 in 1:I){
yi1<-c(yi1,sum(y[u==0 & center==i2]))
yi2<-c(yi2,sum(y[u==1 & center==i2]))
}
##method2##############
fit2 <- glmer(y ~ u + (1|center), family=poisson, data=df,nAGQ=3,verbose =
TRUE, REML = FALSE )
#fit2 <- glmer(y ~ u + (1|center), family=poisson, data=df,nAGQ=3)
estimate_sigma2<-sigma.hat(fit2)$sigma$center
b_estimate<-f...@fixef
alpha1=b_estimate[1]
alpha2=b_estimate[1]+b_estimate[2]
T1=exp(alpha1+estimate_sigma2/2)
T2=exp(alpha2+estimate_sigma2/2)
T3=exp(estimate_sigma2)-1
T4=exp(0.5*estimate_sigma2)
w=exp(s)
eBLP_glimmix_quad[kk]=mean((T4+n*T3*T4*(yi1/n+yi2/n-T1-T2)/(1+T3*n*(T1+T2))-w)^2)
}
MSE_eBLP_quad[i3]<-mean(eBLP_glimmix_quad)
}
print(MSE_eBLP_quad)
then you can execute fit2 again after you run my code. After that, you will
find the following error message.
Error in asMethod(object) : matrix is not symmetric [1,2]
Please let me know how to solve this problem. Thanks a lot.
--
View this message in context:
http://n4.nabble.com/Matrix-is-not-symmetric-under-lme4-tp1836523p1836523.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
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.
