Nothing
R/DPMrandom.R
In DPpackage: Bayesian Nonparametric Modeling in R
### DPMrandom.R
### Extracts random effects from DPpackage objects DPMlmm.
###
### Copyright: Alejandro Jara, 2007-2012,
###
### Last modification: 18-04-2007.
###
### This program is free software; you can redistribute it and/or modify
### it under the terms of the GNU General Public License as published by
### the Free Software Foundation; either version 2 of the License, or (at
### your option) any later version.
###
### This program is distributed in the hope that it will be useful, but
### WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
### General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with this program; if not, write to the Free Software
### Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
###
### The author's contact information:
###
### Alejandro Jara
### Department of Statistics
### Facultad de Matematicas
### Pontificia Universidad Catolica de Chile
### Casilla 306, Correo 22
### Santiago
### Chile
### Voice: +56-2-3544506 URL : http://www.mat.puc.cl/~ajara
### Fax : +56-2-3547729 Email: atjara@uc.cl
###
"DPMrandom"<-
function(object,centered=FALSE,predictive=FALSE,ngrid=1000,gridl=NULL)
UseMethod("DPMrandom")
"DPMrandom.default"<-
function(object,centered=FALSE,predictive=FALSE,ngrid=1000,gridl=NULL)
{
comput<-0
if(is(object, "DPMlmm") || is(object, "DPMolmm") || is(object, "DPMglmm"))
{
comput<-1
random<-matrix(0,nrow=object$nsubject,ncol=object$nrandom)
predp<-rep(0,object$nrandom)
predsd<-rep(0,object$nrandom)
predse<-rep(0,object$nrandom)
predl<-rep(0,object$nrandom)
predu<-rep(0,object$nrandom)
predm<-rep(0,object$nrandom)
randommat<-matrix(object$save.state$randsave,
nrow=object$mcmc$nsave,ncol=object$nrandom*(object$nsubject+1))
dimnames(randommat)<-dimnames(object$save.state$randsave)
thetamat<-matrix(object$save.state$thetasave,nrow=object$mcmc$nsave,
ncol=object$dimen)
counter<-0
for(i in 1:object$nsubject){
for(j in 1:object$nrandom){
counter<-counter+1
if(centered)
{
random[i,j]<-mean(object$save.state$randsave[,counter]-
object$save.state$thetasave[,j])
}
else
{
random[i,j]<-mean(object$save.state$randsave[,counter])
}
}
}
type<-1
dimnames(random)<-list(object$namesre1,object$namesre2)
z<-list(randomm=random,randommat=randommat,thetamat=thetamat,centered=centered,
predictive=predictive,nsubject=object$nsubject,nrandom=object$nrandom,
modelname=object$modelname,call=object$call,type=type,nsave=object$mcmc$nsave)
if(predictive==TRUE)
{
for(i in 1:object$nrandom)
{
counter<-counter+1
if(centered)
{
predp[i]<-mean(object$save.state$randsave[,counter]-
object$save.state$thetasave[,i])
predm[i]<-median(object$save.state$randsave[,counter]-
object$save.state$thetasave[,i])
predsd[i]<-sqrt(var(object$save.state$randsave[,counter]-
object$save.state$thetasave[,i]))
vec<-object$save.state$randsave[,counter]-object$save.state$thetasave[,i]
n<-length(vec)
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
predl[i]<-a$alow[1]
predu[i]<-a$aupp[1]
predse[i]<-predsd[i]/sqrt(n)
}
else
{
predp[i]<-mean(object$save.state$randsave[,counter])
predm[i]<-median(object$save.state$randsave[,counter])
predsd[i]<-sqrt(var(object$save.state$randsave[,counter]))
vec<-object$save.state$randsave[,counter]
n<-length(vec)
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
predl[i]<-a$alow[1]
predu[i]<-a$aupp[1]
predse[i]<-predsd[i]/sqrt(n)
}
}
predtable <- cbind(predp, predm, predsd, predse , predl , predu)
dimnames(predtable) <- list(object$namesre2, c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
z$prediction<-predtable
}
}
if(is(object, "DPMmeta"))
{
comput<-1
if (centered)
{
stop("This option is not implemented for DPMmeta.\n")
}
randommat<-matrix(object$save.state$randsave,
nrow=object$mcmc$nsave,ncol=object$nrec+1)
dimnames(randommat)<-dimnames(object$save.state$randsave)
random<-apply(object$save.state$randsave,2,mean)
random<-random[-(object$nrec+1)]
random<-as.matrix(random,ncol=1)
colnames(random)<-names(object$coefficients)[1]
type<-2
z<-list(randomm=random,modelname=object$modelname,call=object$call,predictive=predictive,
nsubject=object$nrec,nrandom=1,centered=FALSE,randommat=randommat,
type=type,nsave=object$mcmc$nsave)
if(predictive==TRUE)
{
predp<-mean(object$save.state$randsave[,object$nrec+1])
predm<-median(object$save.state$randsave[,object$nrec+1])
predsd<-sqrt(var(object$save.state$randsave[,object$nrec+1]))
vec<-object$save.state$randsave[,object$nrec+1]
n<-length(vec)
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
predl<-a$alow[1]
predu<-a$aupp[1]
predse<-predsd/sqrt(n)
predtable <- cbind(predp, predm, predsd, predse , predl , predu)
dimnames(predtable) <- list("theta", c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
z$prediction<-predtable
}
}
if(comput==1 && predictive==TRUE)
{
q<-object$nrandom
if(q<=2)
{
clustsave<-object$save.state$clustsave
musave<-object$save.state$musave
nsubject<-object$nsubject
nfixed<-object$nfixed
dimen1<-object$nrandom+object$nfixed
dimen2<-0
if(is(object, "DPMlmm"))dimen2<-1
if(is(object, "DPMglmm"))dimen2<-object$dispp
if(is(object, "DPMmeta"))dimen2<-0
dimen3<-object$nrandom*(object$nrandom+1)/2
total<-dimen1+dimen2
sigmakmat<-object$save.state$thetasave[,(total+1):(total+dimen3)]
dimen4<-object$nrandom
total<-dimen1+dimen2+dimen3
mubmat<-object$save.state$thetasave[,(total+1):(total+dimen4)]
dimen5<-object$nrandom*(object$nrandom+1)/2
total<-dimen1+dimen2+dimen3+dimen4
sigmabmat<-object$save.state$thetasave[,(total+1):(total+dimen5)]
total<-dimen1+dimen2+dimen3+dimen4+dimen5
if(is(object, "DPMolmm"))
{
total<-dimen1+dimen2+dimen3+dimen4+dimen5+object$ncateg-2
}
nclusvec<-object$save.state$thetasave[,(total+1)]
alphavec<-object$save.state$thetasave[,(total+2)]
nsave<-length(alphavec)
if(q==1)
{
fs<-rep(0,ngrid)
left<-min(z$randomm)-2*sqrt(var(z$randomm))
right<-max(z$randomm)+2*sqrt(var(z$randomm))
if(is.null(gridl))
{
grid<-seq(left,right,length=ngrid)
}
else
{
gridl<-matrix(gridl,nrow=2)
if(ncol(gridl)!=1)
{
stop("Incorrect number of grid points limits")
}
grid<-seq(gridl[1],gridl[2],length=ngrid)
}
seed1<-sample(1:29000,1)
seed2<-sample(1:29000,1)
foo <- .Fortran("predictivedpmu",
nsave =as.integer(nsave),
nsubject =as.integer(nsubject),
q =as.integer(q),
musave =as.double(musave),
sigmakmat =as.double(sigmakmat),
clustsave =as.integer(clustsave),
nclusvec =as.integer(nclusvec),
alphavec =as.double(alphavec),
mubmat =as.double(mubmat),
sigmabmat =as.double(sigmabmat),
ngrid =as.integer(ngrid),
grid =as.double(grid),
fs =as.double(fs),
seed1 =as.integer(seed1),
seed2 =as.integer(seed2),
PACKAGE ="DPpackage")
z$dens<-foo$fs
z$grid1<-foo$grid
}
if(q==2)
{
ngrid1<-as.integer(sqrt(ngrid))
ngrid2<-ngrid1
f1<-rep(0,ngrid1)
f2<-rep(0,ngrid2)
fs<-matrix(0,nrow=ngrid1,ncol=ngrid2)
if(is.null(gridl))
{
left<-rep(0,2)
right<-rep(0,2)
left[1]<-min(z$randomm[,1])-2*sqrt(var(z$randomm[,1]))
left[2]<-min(z$randomm[,2])-2*sqrt(var(z$randomm[,2]))
right[1]<-max(z$randomm[,1])+2*sqrt(var(z$randomm[,1]))
right[2]<-max(z$randomm[,2])+2*sqrt(var(z$randomm[,2]))
grid1<-seq(left[1],right[1],length=ngrid1)
grid2<-seq(left[2],right[2],length=ngrid1)
}
else
{
gridl<-matrix(gridl,nrow=2)
if(ncol(gridl)!=2)
{
stop("Incorrect number of grid points limits")
}
grid1<-seq(gridl[1,1],gridl[2,1],length=ngrid1)
grid2<-seq(gridl[1,2],gridl[2,2],length=ngrid2)
}
seed1<-sample(1:29000,1)
seed2<-sample(1:29000,1)
iflagr<-rep(0,q)
meanc<-rep(0,q)
meanb<-rep(0,q)
meanw<-rep(0,q)
sigmab<-matrix(0,q,q)
sigmak<-matrix(0,q,q)
theta<-rep(0,q)
workmhr<-rep(0,q*(q+1)/2)
workvr<-rep(0,q)
foo <- .Fortran("predictivedpmb",
nsave =as.integer(nsave),
nsubject =as.integer(nsubject),
q =as.integer(q),
musave =as.double(musave),
sigmakmat =as.double(sigmakmat),
clustsave =as.integer(clustsave),
nclusvec =as.integer(nclusvec),
alphavec =as.double(alphavec),
mubmat =as.double(mubmat),
sigmabmat =as.double(sigmabmat),
ngrid1 =as.integer(ngrid1),
ngrid2 =as.integer(ngrid2),
grid1 =as.double(grid1),
grid2 =as.double(grid2),
fs =as.double(fs),
f1 =as.double(f1),
f2 =as.double(f2),
seed1 =as.integer(seed1),
seed2 =as.integer(seed2),
iflagr =as.integer(iflagr),
meanc =as.double(meanc),
meanb =as.double(meanb),
meanw =as.double(meanw),
sigmab =as.double(sigmab),
sigmak =as.double(sigmak),
theta =as.double(theta),
workmhr =as.double(workmhr),
workvr =as.double(workvr),
PACKAGE ="DPpackage")
z$dens<-matrix(foo$fs,nrow=ngrid1,ncol=ngrid2)
z$grid1<-foo$grid1
z$grid2<-foo$grid2
z$f1<-foo$f1
z$f2<-foo$f2
}
}
}
class(z)<-c("DPMrandom")
return(z)
}
###
### Tools for DPMrandom: print, plot
###
### Copyright: Alejandro Jara, 2007
### Last modification: 02-04-2007.
"print.DPMrandom"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n","Random effect information for the DP object:","\n\nCall:\n",sep = "")
print(x$call)
cat("\n")
if(x$predictive)
{
cat("\nPredictive distribution:\n")
print.default(format(x$prediction, digits = digits), print.gap = 2,
quote = FALSE)
}
else
{
cat("\nPosterior mean of subject-specific components:\n\n")
print.default(format(x$randomm, digits = digits), print.gap = 2,
quote = FALSE)
}
cat("\n\n")
invisible(x)
}
"plot.DPMrandom"<-function(x, ask=TRUE, hpd=TRUE, nfigr=2, nfigc=2, subject=NULL, col="#bdfcc9", ...)
{
fancydensplot<-function(x, hpd=TRUE, npts=200,
xlim=NULL,xlab="", ylab="", main="",col="#bdfcc9", ...)
# Author: AJV, 2006
#
{
dens <- density(x,n=npts)
densx <- dens$x
densy <- dens$y
meanvar <- mean(x)
densx1 <- max(densx[densx<=meanvar])
densx2 <- min(densx[densx>=meanvar])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ymean <- densy1 + ((densy2-densy1)/(densx2-densx1))*(meanvar-densx1)
if(hpd==TRUE)
{
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
n<-length(x)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(x),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
xlinf<-a$alow[1]
xlsup<-a$aupp[1]
}
else
{
xlinf <- quantile(x,0.025)
xlsup <- quantile(x,0.975)
}
densx1 <- max(densx[densx<=xlinf])
densx2 <- min(densx[densx>=xlinf])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylinf <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlinf-densx1)
densx1 <- max(densx[densx<=xlsup])
densx2 <- min(densx[densx>=xlsup])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylsup <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlsup-densx1)
if(is.null(xlim))
{
xlim<-c(min(densx), max(densx))
}
plot(0.,0.,xlim = xlim, ylim = c(min(densy), max(densy)),
axes = F,type = "n" , xlab=xlab, ylab=ylab, main=main, cex=1.2)
xpol<-c(xlinf,xlinf,densx[densx>=xlinf & densx <=xlsup],xlsup,xlsup)
ypol<-c(0,ylinf,densy[densx>=xlinf & densx <=xlsup] ,ylsup,0)
polygon(xpol, ypol, border = FALSE,col=col)
lines(xlim,c(0,0),lwd=1.2)
segments(xlim[1],0, xlim[1],max(densy),lwd=1.2)
lines(densx,densy,lwd=1.2)
segments(meanvar, 0, meanvar, ymean,lwd=1.2)
segments(xlinf, 0, xlinf, ylinf,lwd=1.2)
segments(xlsup, 0, xlsup, ylsup,lwd=1.2)
axis(1., at = round(c(xlinf, meanvar,xlsup), 2.), labels = T,pos = 0.)
axis(1., at = round(seq(xlim[1],xlim[2],length=15), 2.), labels = F,pos = 0.)
axis(2., at = round(seq(0,max(densy),length=5), 2.), labels = T,pos =xlim[1])
}
if(is(x, "DPMrandom")){
oldpar <- par(no.readonly = TRUE)
par(ask = ask)
if(x$predictive)
{
coef.p<-x$randomm
n<-dim(coef.p)[2]
pnames<-colnames(coef.p)
start<-x$nsubject*n
if(n==1)
{
nfigr<-1
nfigc<-1
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
xx<-matrix(x$grid1,ncol=1)
dens<-x$dens
title<-paste("Density of",pnames[1],sep=" ")
plot(0.,0.,xlim = c(min(xx),max(xx)), ylim = c(0, (max(dens)+0.01)),
axes = F,type = "n" , xlab="values", ylab="density", main=title, cex=1.2)
lines(x$grid1,x$dens,lwd=2)
lines(c(min(xx),max(xx)),c(0,0),lwd=1.2)
segments(min(xx),0, min(xx),max(dens)+0.01,lwd=1.2)
axis(1., at = round(seq(min(xx),max(xx),length=15), 2.), labels = T,pos = 0.)
axis(2., at = round(seq(0,max(dens)+0.01,length=5), 2.), labels = T,pos =min(xx))
for(j in 1:x$nsubject)
{
points(x$randomm[j,1],0,col="red",pch=20)
}
}
if(n==2)
{
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
xx<-matrix(x$grid1,ncol=1)
yy<-matrix(x$grid2,ncol=1)
z<-x$dens
ngrid<-length(xx)
dens1<-x$f1
dens2<-x$f2
colnames(xx)<-pnames[1]
colnames(yy)<-pnames[2]
title<-paste("Density of",pnames[1],sep=" ")
plot(0.,0.,xlim = c(min(xx),max(xx)), ylim = c(0, (max(dens1)+0.01)),
axes = F,type = "n" , xlab="values", ylab="density", main=title, cex=1.2)
lines(xx,dens1,lwd=2)
lines(c(min(xx),max(xx)),c(0,0),lwd=1.2)
segments(min(xx),0, min(xx),max(dens1)+0.01,lwd=1.2)
axis(1., at = round(seq(min(xx),max(xx),length=15), 2.), labels = T,pos = 0.)
axis(2., at = round(seq(0,max(dens1)+0.01,length=5), 2.), labels = T,pos =min(xx))
for(j in 1:x$nsubject)
{
points(x$randomm[j,1],0,col="red",pch=20)
}
title<-paste("Density of",pnames[2],sep=" ")
plot(0.,0.,xlim = c(min(yy),max(yy)), ylim = c(0, (max(dens2)+0.01)),
axes = F,type = "n" , xlab="values", ylab="density", main=title, cex=1.2)
lines(yy,dens2,lwd=2)
lines(c(min(yy),max(yy)),c(0,0),lwd=1.2)
segments(min(yy),0, min(yy),max(dens2)+0.01,lwd=1.2)
axis(1., at = round(seq(min(yy),max(yy),length=15), 2.), labels = T,pos = 0.)
axis(2., at = round(seq(0,max(dens2)+0.01,length=5), 2.), labels = T,pos =min(yy))
for(j in 1:x$nsubject)
{
points(x$randomm[j,2],0,col="red",pch=20)
}
title<-paste("Density of",pnames[1],sep=" ")
title<-paste(title,pnames[2],sep="-")
contour(xx,yy,z,xlab=pnames[1],ylab=pnames[2],main=title)
points(x$randomm[,1],x$randomm[,2],col="black",pch=20)
persp(xx,yy,z,xlab=pnames[1],ylab=pnames[2],zlab="density",theta=-30,phi=15,expand = 0.9, ltheta = 120,main=title)
}
if(n>2)
{
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
for(i in 1:n)
{
title1<-paste("Trace of",pnames[i],sep=" ")
title2<-paste("Density of",pnames[i],sep=" ")
if(x$centered)
{
vec<-(x$randommat[,start+i]-x$thetamat[,i])
}
else vec<-x$randommat[,start+i]
names(vec)<-pnames[i]
vectmp<-x$randomm[,i]
xlim<-c(min(vectmp)-6.5*sqrt(var(vectmp)),max(vectmp)+6.5*sqrt(var(vectmp)))
plot(vec,type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(vec,xlim=xlim,hpd=hpd,main=title2,xlab="values", ylab="density", col=col)
for(j in 1:x$nsubject)
{
points(x$randomm[j,i],0,col="red",pch=20)
}
}
}
}
else
{
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
coef.p<-x$randommat
n<-dim(coef.p)[2]
pnames<-colnames(coef.p)
count<-0
for(i in 1:x$nsubject)
{
for(j in 1:x$nrandom)
{
count<-count+1
title1<-paste("Trace of",pnames[count],sep=" ")
title2<-paste("Density of",pnames[count],sep=" ")
if(x$centered)
{
vec<-(x$randommat[,count]-x$thetamat[,j])
}
else vec<-x$randommat[,count]
plot(vec,type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(vec,hpd=hpd,main=title2,xlab="values", ylab="density", col=col)
}
}
}
par(oldpar)
}
}
Try the DPpackage package in your browser
Any scripts or data that you put into this service are public.
DPpackage documentation built on May 1, 2019, 10:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### DPMrandom.R
### Extracts random effects from DPpackage objects DPMlmm.
###
### Copyright: Alejandro Jara, 2007-2012,
###
### Last modification: 18-04-2007.
###
### This program is free software; you can redistribute it and/or modify
### it under the terms of the GNU General Public License as published by
### the Free Software Foundation; either version 2 of the License, or (at
### your option) any later version.
###
### This program is distributed in the hope that it will be useful, but
### WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
### General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with this program; if not, write to the Free Software
### Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
###
### The author's contact information:
###
### Alejandro Jara
### Department of Statistics
### Facultad de Matematicas
### Pontificia Universidad Catolica de Chile
### Casilla 306, Correo 22
### Santiago
### Chile
### Voice: +56-2-3544506 URL : http://www.mat.puc.cl/~ajara
### Fax : +56-2-3547729 Email: atjara@uc.cl
###
"DPMrandom"<-
function(object,centered=FALSE,predictive=FALSE,ngrid=1000,gridl=NULL)
UseMethod("DPMrandom")
"DPMrandom.default"<-
function(object,centered=FALSE,predictive=FALSE,ngrid=1000,gridl=NULL)
{
comput<-0
if(is(object, "DPMlmm") || is(object, "DPMolmm") || is(object, "DPMglmm"))
{
comput<-1
random<-matrix(0,nrow=object$nsubject,ncol=object$nrandom)
predp<-rep(0,object$nrandom)
predsd<-rep(0,object$nrandom)
predse<-rep(0,object$nrandom)
predl<-rep(0,object$nrandom)
predu<-rep(0,object$nrandom)
predm<-rep(0,object$nrandom)
randommat<-matrix(object$save.state$randsave,
nrow=object$mcmc$nsave,ncol=object$nrandom*(object$nsubject+1))
dimnames(randommat)<-dimnames(object$save.state$randsave)
thetamat<-matrix(object$save.state$thetasave,nrow=object$mcmc$nsave,
ncol=object$dimen)
counter<-0
for(i in 1:object$nsubject){
for(j in 1:object$nrandom){
counter<-counter+1
if(centered)
{
random[i,j]<-mean(object$save.state$randsave[,counter]-
object$save.state$thetasave[,j])
}
else
{
random[i,j]<-mean(object$save.state$randsave[,counter])
}
}
}
type<-1
dimnames(random)<-list(object$namesre1,object$namesre2)
z<-list(randomm=random,randommat=randommat,thetamat=thetamat,centered=centered,
predictive=predictive,nsubject=object$nsubject,nrandom=object$nrandom,
modelname=object$modelname,call=object$call,type=type,nsave=object$mcmc$nsave)
if(predictive==TRUE)
{
for(i in 1:object$nrandom)
{
counter<-counter+1
if(centered)
{
predp[i]<-mean(object$save.state$randsave[,counter]-
object$save.state$thetasave[,i])
predm[i]<-median(object$save.state$randsave[,counter]-
object$save.state$thetasave[,i])
predsd[i]<-sqrt(var(object$save.state$randsave[,counter]-
object$save.state$thetasave[,i]))
vec<-object$save.state$randsave[,counter]-object$save.state$thetasave[,i]
n<-length(vec)
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
predl[i]<-a$alow[1]
predu[i]<-a$aupp[1]
predse[i]<-predsd[i]/sqrt(n)
}
else
{
predp[i]<-mean(object$save.state$randsave[,counter])
predm[i]<-median(object$save.state$randsave[,counter])
predsd[i]<-sqrt(var(object$save.state$randsave[,counter]))
vec<-object$save.state$randsave[,counter]
n<-length(vec)
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
predl[i]<-a$alow[1]
predu[i]<-a$aupp[1]
predse[i]<-predsd[i]/sqrt(n)
}
}
predtable <- cbind(predp, predm, predsd, predse , predl , predu)
dimnames(predtable) <- list(object$namesre2, c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
z$prediction<-predtable
}
}
if(is(object, "DPMmeta"))
{
comput<-1
if (centered)
{
stop("This option is not implemented for DPMmeta.\n")
}
randommat<-matrix(object$save.state$randsave,
nrow=object$mcmc$nsave,ncol=object$nrec+1)
dimnames(randommat)<-dimnames(object$save.state$randsave)
random<-apply(object$save.state$randsave,2,mean)
random<-random[-(object$nrec+1)]
random<-as.matrix(random,ncol=1)
colnames(random)<-names(object$coefficients)[1]
type<-2
z<-list(randomm=random,modelname=object$modelname,call=object$call,predictive=predictive,
nsubject=object$nrec,nrandom=1,centered=FALSE,randommat=randommat,
type=type,nsave=object$mcmc$nsave)
if(predictive==TRUE)
{
predp<-mean(object$save.state$randsave[,object$nrec+1])
predm<-median(object$save.state$randsave[,object$nrec+1])
predsd<-sqrt(var(object$save.state$randsave[,object$nrec+1]))
vec<-object$save.state$randsave[,object$nrec+1]
n<-length(vec)
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
predl<-a$alow[1]
predu<-a$aupp[1]
predse<-predsd/sqrt(n)
predtable <- cbind(predp, predm, predsd, predse , predl , predu)
dimnames(predtable) <- list("theta", c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
z$prediction<-predtable
}
}
if(comput==1 && predictive==TRUE)
{
q<-object$nrandom
if(q<=2)
{
clustsave<-object$save.state$clustsave
musave<-object$save.state$musave
nsubject<-object$nsubject
nfixed<-object$nfixed
dimen1<-object$nrandom+object$nfixed
dimen2<-0
if(is(object, "DPMlmm"))dimen2<-1
if(is(object, "DPMglmm"))dimen2<-object$dispp
if(is(object, "DPMmeta"))dimen2<-0
dimen3<-object$nrandom*(object$nrandom+1)/2
total<-dimen1+dimen2
sigmakmat<-object$save.state$thetasave[,(total+1):(total+dimen3)]
dimen4<-object$nrandom
total<-dimen1+dimen2+dimen3
mubmat<-object$save.state$thetasave[,(total+1):(total+dimen4)]
dimen5<-object$nrandom*(object$nrandom+1)/2
total<-dimen1+dimen2+dimen3+dimen4
sigmabmat<-object$save.state$thetasave[,(total+1):(total+dimen5)]
total<-dimen1+dimen2+dimen3+dimen4+dimen5
if(is(object, "DPMolmm"))
{
total<-dimen1+dimen2+dimen3+dimen4+dimen5+object$ncateg-2
}
nclusvec<-object$save.state$thetasave[,(total+1)]
alphavec<-object$save.state$thetasave[,(total+2)]
nsave<-length(alphavec)
if(q==1)
{
fs<-rep(0,ngrid)
left<-min(z$randomm)-2*sqrt(var(z$randomm))
right<-max(z$randomm)+2*sqrt(var(z$randomm))
if(is.null(gridl))
{
grid<-seq(left,right,length=ngrid)
}
else
{
gridl<-matrix(gridl,nrow=2)
if(ncol(gridl)!=1)
{
stop("Incorrect number of grid points limits")
}
grid<-seq(gridl[1],gridl[2],length=ngrid)
}
seed1<-sample(1:29000,1)
seed2<-sample(1:29000,1)
foo <- .Fortran("predictivedpmu",
nsave =as.integer(nsave),
nsubject =as.integer(nsubject),
q =as.integer(q),
musave =as.double(musave),
sigmakmat =as.double(sigmakmat),
clustsave =as.integer(clustsave),
nclusvec =as.integer(nclusvec),
alphavec =as.double(alphavec),
mubmat =as.double(mubmat),
sigmabmat =as.double(sigmabmat),
ngrid =as.integer(ngrid),
grid =as.double(grid),
fs =as.double(fs),
seed1 =as.integer(seed1),
seed2 =as.integer(seed2),
PACKAGE ="DPpackage")
z$dens<-foo$fs
z$grid1<-foo$grid
}
if(q==2)
{
ngrid1<-as.integer(sqrt(ngrid))
ngrid2<-ngrid1
f1<-rep(0,ngrid1)
f2<-rep(0,ngrid2)
fs<-matrix(0,nrow=ngrid1,ncol=ngrid2)
if(is.null(gridl))
{
left<-rep(0,2)
right<-rep(0,2)
left[1]<-min(z$randomm[,1])-2*sqrt(var(z$randomm[,1]))
left[2]<-min(z$randomm[,2])-2*sqrt(var(z$randomm[,2]))
right[1]<-max(z$randomm[,1])+2*sqrt(var(z$randomm[,1]))
right[2]<-max(z$randomm[,2])+2*sqrt(var(z$randomm[,2]))
grid1<-seq(left[1],right[1],length=ngrid1)
grid2<-seq(left[2],right[2],length=ngrid1)
}
else
{
gridl<-matrix(gridl,nrow=2)
if(ncol(gridl)!=2)
{
stop("Incorrect number of grid points limits")
}
grid1<-seq(gridl[1,1],gridl[2,1],length=ngrid1)
grid2<-seq(gridl[1,2],gridl[2,2],length=ngrid2)
}
seed1<-sample(1:29000,1)
seed2<-sample(1:29000,1)
iflagr<-rep(0,q)
meanc<-rep(0,q)
meanb<-rep(0,q)
meanw<-rep(0,q)
sigmab<-matrix(0,q,q)
sigmak<-matrix(0,q,q)
theta<-rep(0,q)
workmhr<-rep(0,q*(q+1)/2)
workvr<-rep(0,q)
foo <- .Fortran("predictivedpmb",
nsave =as.integer(nsave),
nsubject =as.integer(nsubject),
q =as.integer(q),
musave =as.double(musave),
sigmakmat =as.double(sigmakmat),
clustsave =as.integer(clustsave),
nclusvec =as.integer(nclusvec),
alphavec =as.double(alphavec),
mubmat =as.double(mubmat),
sigmabmat =as.double(sigmabmat),
ngrid1 =as.integer(ngrid1),
ngrid2 =as.integer(ngrid2),
grid1 =as.double(grid1),
grid2 =as.double(grid2),
fs =as.double(fs),
f1 =as.double(f1),
f2 =as.double(f2),
seed1 =as.integer(seed1),
seed2 =as.integer(seed2),
iflagr =as.integer(iflagr),
meanc =as.double(meanc),
meanb =as.double(meanb),
meanw =as.double(meanw),
sigmab =as.double(sigmab),
sigmak =as.double(sigmak),
theta =as.double(theta),
workmhr =as.double(workmhr),
workvr =as.double(workvr),
PACKAGE ="DPpackage")
z$dens<-matrix(foo$fs,nrow=ngrid1,ncol=ngrid2)
z$grid1<-foo$grid1
z$grid2<-foo$grid2
z$f1<-foo$f1
z$f2<-foo$f2
}
}
}
class(z)<-c("DPMrandom")
return(z)
}
###
### Tools for DPMrandom: print, plot
###
### Copyright: Alejandro Jara, 2007
### Last modification: 02-04-2007.
"print.DPMrandom"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n","Random effect information for the DP object:","\n\nCall:\n",sep = "")
print(x$call)
cat("\n")
if(x$predictive)
{
cat("\nPredictive distribution:\n")
print.default(format(x$prediction, digits = digits), print.gap = 2,
quote = FALSE)
}
else
{
cat("\nPosterior mean of subject-specific components:\n\n")
print.default(format(x$randomm, digits = digits), print.gap = 2,
quote = FALSE)
}
cat("\n\n")
invisible(x)
}
"plot.DPMrandom"<-function(x, ask=TRUE, hpd=TRUE, nfigr=2, nfigc=2, subject=NULL, col="#bdfcc9", ...)
{
fancydensplot<-function(x, hpd=TRUE, npts=200,
xlim=NULL,xlab="", ylab="", main="",col="#bdfcc9", ...)
# Author: AJV, 2006
#
{
dens <- density(x,n=npts)
densx <- dens$x
densy <- dens$y
meanvar <- mean(x)
densx1 <- max(densx[densx<=meanvar])
densx2 <- min(densx[densx>=meanvar])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ymean <- densy1 + ((densy2-densy1)/(densx2-densx1))*(meanvar-densx1)
if(hpd==TRUE)
{
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
n<-length(x)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(x),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
xlinf<-a$alow[1]
xlsup<-a$aupp[1]
}
else
{
xlinf <- quantile(x,0.025)
xlsup <- quantile(x,0.975)
}
densx1 <- max(densx[densx<=xlinf])
densx2 <- min(densx[densx>=xlinf])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylinf <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlinf-densx1)
densx1 <- max(densx[densx<=xlsup])
densx2 <- min(densx[densx>=xlsup])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylsup <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlsup-densx1)
if(is.null(xlim))
{
xlim<-c(min(densx), max(densx))
}
plot(0.,0.,xlim = xlim, ylim = c(min(densy), max(densy)),
axes = F,type = "n" , xlab=xlab, ylab=ylab, main=main, cex=1.2)
xpol<-c(xlinf,xlinf,densx[densx>=xlinf & densx <=xlsup],xlsup,xlsup)
ypol<-c(0,ylinf,densy[densx>=xlinf & densx <=xlsup] ,ylsup,0)
polygon(xpol, ypol, border = FALSE,col=col)
lines(xlim,c(0,0),lwd=1.2)
segments(xlim[1],0, xlim[1],max(densy),lwd=1.2)
lines(densx,densy,lwd=1.2)
segments(meanvar, 0, meanvar, ymean,lwd=1.2)
segments(xlinf, 0, xlinf, ylinf,lwd=1.2)
segments(xlsup, 0, xlsup, ylsup,lwd=1.2)
axis(1., at = round(c(xlinf, meanvar,xlsup), 2.), labels = T,pos = 0.)
axis(1., at = round(seq(xlim[1],xlim[2],length=15), 2.), labels = F,pos = 0.)
axis(2., at = round(seq(0,max(densy),length=5), 2.), labels = T,pos =xlim[1])
}
if(is(x, "DPMrandom")){
oldpar <- par(no.readonly = TRUE)
par(ask = ask)
if(x$predictive)
{
coef.p<-x$randomm
n<-dim(coef.p)[2]
pnames<-colnames(coef.p)
start<-x$nsubject*n
if(n==1)
{
nfigr<-1
nfigc<-1
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
xx<-matrix(x$grid1,ncol=1)
dens<-x$dens
title<-paste("Density of",pnames[1],sep=" ")
plot(0.,0.,xlim = c(min(xx),max(xx)), ylim = c(0, (max(dens)+0.01)),
axes = F,type = "n" , xlab="values", ylab="density", main=title, cex=1.2)
lines(x$grid1,x$dens,lwd=2)
lines(c(min(xx),max(xx)),c(0,0),lwd=1.2)
segments(min(xx),0, min(xx),max(dens)+0.01,lwd=1.2)
axis(1., at = round(seq(min(xx),max(xx),length=15), 2.), labels = T,pos = 0.)
axis(2., at = round(seq(0,max(dens)+0.01,length=5), 2.), labels = T,pos =min(xx))
for(j in 1:x$nsubject)
{
points(x$randomm[j,1],0,col="red",pch=20)
}
}
if(n==2)
{
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
xx<-matrix(x$grid1,ncol=1)
yy<-matrix(x$grid2,ncol=1)
z<-x$dens
ngrid<-length(xx)
dens1<-x$f1
dens2<-x$f2
colnames(xx)<-pnames[1]
colnames(yy)<-pnames[2]
title<-paste("Density of",pnames[1],sep=" ")
plot(0.,0.,xlim = c(min(xx),max(xx)), ylim = c(0, (max(dens1)+0.01)),
axes = F,type = "n" , xlab="values", ylab="density", main=title, cex=1.2)
lines(xx,dens1,lwd=2)
lines(c(min(xx),max(xx)),c(0,0),lwd=1.2)
segments(min(xx),0, min(xx),max(dens1)+0.01,lwd=1.2)
axis(1., at = round(seq(min(xx),max(xx),length=15), 2.), labels = T,pos = 0.)
axis(2., at = round(seq(0,max(dens1)+0.01,length=5), 2.), labels = T,pos =min(xx))
for(j in 1:x$nsubject)
{
points(x$randomm[j,1],0,col="red",pch=20)
}
title<-paste("Density of",pnames[2],sep=" ")
plot(0.,0.,xlim = c(min(yy),max(yy)), ylim = c(0, (max(dens2)+0.01)),
axes = F,type = "n" , xlab="values", ylab="density", main=title, cex=1.2)
lines(yy,dens2,lwd=2)
lines(c(min(yy),max(yy)),c(0,0),lwd=1.2)
segments(min(yy),0, min(yy),max(dens2)+0.01,lwd=1.2)
axis(1., at = round(seq(min(yy),max(yy),length=15), 2.), labels = T,pos = 0.)
axis(2., at = round(seq(0,max(dens2)+0.01,length=5), 2.), labels = T,pos =min(yy))
for(j in 1:x$nsubject)
{
points(x$randomm[j,2],0,col="red",pch=20)
}
title<-paste("Density of",pnames[1],sep=" ")
title<-paste(title,pnames[2],sep="-")
contour(xx,yy,z,xlab=pnames[1],ylab=pnames[2],main=title)
points(x$randomm[,1],x$randomm[,2],col="black",pch=20)
persp(xx,yy,z,xlab=pnames[1],ylab=pnames[2],zlab="density",theta=-30,phi=15,expand = 0.9, ltheta = 120,main=title)
}
if(n>2)
{
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
for(i in 1:n)
{
title1<-paste("Trace of",pnames[i],sep=" ")
title2<-paste("Density of",pnames[i],sep=" ")
if(x$centered)
{
vec<-(x$randommat[,start+i]-x$thetamat[,i])
}
else vec<-x$randommat[,start+i]
names(vec)<-pnames[i]
vectmp<-x$randomm[,i]
xlim<-c(min(vectmp)-6.5*sqrt(var(vectmp)),max(vectmp)+6.5*sqrt(var(vectmp)))
plot(vec,type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(vec,xlim=xlim,hpd=hpd,main=title2,xlab="values", ylab="density", col=col)
for(j in 1:x$nsubject)
{
points(x$randomm[j,i],0,col="red",pch=20)
}
}
}
}
else
{
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr , ncol=nfigc ,byrow=TRUE))
coef.p<-x$randommat
n<-dim(coef.p)[2]
pnames<-colnames(coef.p)
count<-0
for(i in 1:x$nsubject)
{
for(j in 1:x$nrandom)
{
count<-count+1
title1<-paste("Trace of",pnames[count],sep=" ")
title2<-paste("Density of",pnames[count],sep=" ")
if(x$centered)
{
vec<-(x$randommat[,count]-x$thetamat[,j])
}
else vec<-x$randommat[,count]
plot(vec,type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(vec,hpd=hpd,main=title2,xlab="values", ylab="density", col=col)
}
}
}
par(oldpar)
}
}
Try the DPpackage package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.