# NOT RUN {
##############################################
# A simulated data using "perfect"
# simulation from a mixture of two
# normals and normal true models for
# the random effects.
# A Poisson sampling distribution
# is considered.
##############################################
# Functions needed to simulate random effects
# and to evaluate true models
findq <- function(true.cdf,target,low,
upp,epsilon=0.0000001)
{
plow <- true.cdf(low)
pupp <- true.cdf(upp)
pcenter <- true.cdf((upp+low)/2)
err <- abs(pcenter-target)
i <- 0
while(err > epsilon)
{
i <- i + 1
if(target< pcenter)
{
upp <- (upp+low)/2
pupp <- pcenter
pcenter <- true.cdf((upp+low)/2)
err <- abs(pcenter-target)
}
if(target>= pcenter)
{
low <- (upp+low)/2
plow <- pcenter
pcenter <- true.cdf((upp+low)/2)
err <- abs(pcenter-target)
}
}
return((upp+low)/2)
}
true.dens1 <- function(x)
{
0.5*dnorm(x,2.,sqrt(0.005))+
0.5*dnorm(x,2.85,sqrt(0.005))
}
true.dens2 <- function(x)
{
dnorm(x,2.1,sqrt(0.0324))
}
true.cdf1 <- function(x)
{
0.5*pnorm(x,2.,sqrt(0.005))+
0.5*pnorm(x,2.85,sqrt(0.005))
}
true.cdf2 <- function(x)
{
pnorm(x,2.1,sqrt(0.0324))
}
# Simulation of random effects
nsubject <- 200
nsim <- nsubject/2
qq <- seq(1,nsim)/(nsim+1)
b1 <- rep(0,nsim)
for(i in 1:nsim)
{
aa <- findq(true.cdf1,qq[i],low=-6,upp=6)
b1[i] <- aa
}
b2 <- rep(0,nsim)
for(i in 1:nsim)
{
aa <- findq(true.cdf2,qq[i],low=-6,upp=6)
b2[i] <- aa
}
trt <- c(rep(0,nsim),rep(1,nsim))
b <- c(b1,b2)
xtf <- cbind(rep(1,nsubject),trt)
# Simulation of responses
ni <- 5
nrec <- nsubject*ni
y <- NULL
g <- NULL
x <- NULL
z <- rnorm(nrec)
ll <- 0
for(i in 1:nsubject)
{
g <- c(g,rep(i,ni))
for(j in 1:ni)
{
ll <- ll +1
etaij <- b[i] + 1.2*z[ll]
ytmp <- rpois(1,exp(etaij))
y <- c(y,ytmp)
x <- rbind(x,c(1,trt[i],z[ll]))
}
}
colnames(x) <- c("Intercept","trt","z")
# Design matrix for prediction
xpred <- rbind(c(1,0,0),c(1,1,0))
xtfpred <- rbind(c(1,0),c(1,1))
prediction <- list(xpred=xpred,
xtfpred=xtfpred,
quans=c(0.03,0.50,0.97))
# Prior information
prior <- list(maxm=5,
alpha=0.5,
mub=rep(0,3),
Sb=diag(1000,3),
taub1=2.002,
taub2=2.002)
# Initial state
state <- NULL
# MCMC parameters
nburn <- 4000
nsave <- 4000
nskip <- 3
ndisplay <- 500
mcmc <- list(nburn=nburn,
nsave=nsave,
nskip=nskip,
ndisplay=ndisplay)
# Fitting the model
fit1 <- LDTFPglmm(y=y,x=x,g=g,family=poisson(log),
xtf=xtf,grid=seq(1.2,3.2,len=200),
prediction=prediction,
prior=prior,
mcmc=mcmc,
state=state,
status=TRUE,
compute.band=TRUE)
# Plotting density estimates and true models
# for the random intercepts
par(cex=1.7,mar=c(4.1, 4.1, 1, 1))
plot(fit1$grid,fit1$densu[1,],type="l",xlab="y",
ylab="f(y|x)",lty=2,lwd=3,main="trt=0")
lines(fit1$grid,fit1$densl[1,],lty=2,lwd=3)
lines(fit1$grid,fit1$densm[1,],lty=1,lwd=3)
tmp1 <- true.dens1(fit1$grid)
lines(fit1$grid,tmp1,lty=1,lwd=3,col="red")
par(cex=1.7,mar=c(4.1, 4.1, 1, 1))
plot(fit1$grid,fit1$densu[2,],type="l",xlab="y",
ylab="f(y|x)",lty=2,lwd=3,main="trt=1")
lines(fit1$grid,fit1$densl[2,],lty=2,lwd=3)
lines(fit1$grid,fit1$densm[2,],lty=1,lwd=3)
tmp1 <- true.dens2(fit1$grid)
lines(fit1$grid,tmp1,lty=1,lwd=3,col="red")
# }
Run the code above in your browser using DataLab