# NOT RUN {
########################
# IgG data
########################
data(igg)
z <- igg$age
y <- log(igg$igg)
# Design matrices
ages1 <- z^2
ages2 <- 1/ages1
x <- cbind(rep(1,length(y)),ages1,ages2)
xtf <- cbind(rep(1,length(y)),ages1,ages2)
colnames(x) <- c("(Intercept)","age^2","age^{-2}")
colnames(xtf) <- c("(Intercept)","age^2","age^{-2}")
# Prediction
xdpred <- c(11/12,25/12,38/12,52/12,65/12)
agesp1 <- xdpred^2
agesp2 <- 1/agesp1
xdenpred <- cbind(rep(1,length(xdpred)),agesp1,agesp2)
xtfdenpred <- xdenpred
xmpred <- seq(0.5,6,len=50)
agesp1 <- xmpred^2
agesp2 <- 1/agesp1
xmedpred <- cbind(rep(1,length(xmpred)),agesp1,agesp2)
xtfmedpred <- xmedpred
prediction <- list(xdenpred=xdenpred,
xtfdenpred=xtfdenpred,
xmedpred=xmedpred,
xtfmedpred=xtfmedpred,
quans=c(0.03,0.50,0.97))
# Prior information
Sb <- diag(1000,3)
mub <- rep(0,3)
prior <- list(maxm=4,
a0=1,
b0=1,
mub=mub,
Sb=Sb,
tau1=2.02,
tau2=2.02)
# Initial state
state <- NULL
# MCMC parameters
mcmc <- list(nburn=5000,
nsave=5000,
nskip=4,
ndisplay=200)
# Fitting the model
fit1 <- LDTFPdensity(y=y,
x=x,
xtf=xtf,
prediction=prediction,
prior=prior,
mcmc=mcmc,
state=state,
status=TRUE,
compute.band=TRUE)
fit1
summary(fit1)
plot(fit1)
# Plots predictions
# (conditional densities and quantile functions)
par(mfrow=c(3,2))
for(i in 1:5)
{
plot(fit1$grid,fit1$densm[i,],lwd=2,
type="l",lty=1,col="black",xlab="log IgG",ylab="density",
ylim=c(0,2))
lines(fit1$grid,fit1$densl[i,],lwd=1,lty=2,col="black")
lines(fit1$grid,fit1$densu[i,],lwd=1,lty=2,col="black")
}
plot(z,y,ylab="log IgG",xlab="Age (years)")
lines(xmpred,fit1$qmm[,2],lwd=2)
lines(xmpred,fit1$qml[,2],lwd=1,lty=2)
lines(xmpred,fit1$qmu[,2],lwd=1,lty=2)
lines(xmpred,fit1$qmm[,1],lwd=2)
lines(xmpred,fit1$qml[,1],lwd=1,lty=2)
lines(xmpred,fit1$qmu[,1],lwd=1,lty=2)
lines(xmpred,fit1$qmm[,3],lwd=2)
lines(xmpred,fit1$qml[,3],lwd=1,lty=2)
lines(xmpred,fit1$qmu[,3],lwd=1,lty=2)
#######################################
# A simulated data using "perfect"
# simulation (mixture of two normals
# and normal true models).
#######################################
# Functions needed to simulate data
# 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.5,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.50,sqrt(0.005))+
0.5*pnorm(x,2.85,sqrt(0.005))
}
true.cdf2 <- function(x)
{
pnorm(x,2.1,sqrt(0.0324))
}
# Simulation
nsim <- 500
qq <- seq(1,nsim)/(nsim+1)
y1 <- rep(0,nsim)
for(i in 1:nsim)
{
aa <- findq(true.cdf1,qq[i],low=-6,upp=6)
y1[i] <- aa
}
y2 <- rep(0,nsim)
for(i in 1:nsim)
{
aa <- findq(true.cdf2,qq[i],low=-6,upp=6)
y2[i] <- aa
}
trt <- c(rep(0,nsim),rep(1,nsim))
y <- c(y1,y2)
# Design matrices
W1 <- cbind(rep(1,2*nsim),trt)
W2 <- cbind(rep(1,2*nsim),trt)
colnames(W1) <- c("(Intercept)","trt")
colnames(W2) <- c("(Intercept)","trt")
# Design matrix for prediction
intp <- rep(1,2)
trtp <- c(0,1)
zpred <- cbind(intp,trtp)
prediction <- list(xdenpred=zpred,
xtfdenpred=zpred,
xmedpred=zpred,
xtfmedpred=zpred,
quans=c(0.03,0.50,0.97))
# Prior information
prior <- list(maxm=5,
a0=1,
b0=1,
mub=rep(0,2),
Sb=diag(1000,2),
tau1=2.002,
tau2=2.002)
# Initial state
state <- NULL
# MCMC parameters
nburn <- 5000
nsave <- 5000
nskip <- 4
ndisplay <- 200
mcmc <- list(nburn=nburn,
nsave=nsave,
nskip=nskip,
ndisplay=ndisplay)
# Fitting the model
fit1 <- LDTFPdensity(y=y,
x=W1,
xtf=W2,
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
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