## 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")
#
# ## End(Not run)
Run the code above in your browser using DataLab