https://github.com/cran/HLSM
Tip revision: dcd2fa6f1994b6500075e0b6f0310a45f112b45f authored by Samrachana Adhikari on 18 June 2014, 00:00:00 UTC
version 0.2
version 0.2
Tip revision: dcd2fa6
PlottingHLSMFitsfromHLSMPackage.r
###Plots####
####DENSITY PLOTS######
###function to draw density plots and color the 95% equal tailed credible interval###
#if(FALSE){
#plot.posterior.density=function(chain, title, color){
# plot(density(chain), main=title, ylab="", xlab="")
# Qs=quantile(chain, probs=c(0.025, 0.975, 0.5))
# x1 <- min(which(density(chain)$x >= Qs[1]))
#x2 <- max(which(density(chain)$x < Qs[2]))
#with(density(chain), polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), density=NA, col=color))
#abline(v=Qs[3], col="gray25")
#}
#}
######"BOX" PLOTS########
###chain.list is a list MCMC output, either a vector for one parameter or matrix for several parameters###
plotHLSM.random.fit=function(fitted.model, parameter, burnin = 0, thin = 1){
colors=c("navy", "orange", "darkmagenta")
# niter = fitted.model$call$niter
# if(is.null(burnin)){burnin = round(0.1*niter)}
# if(is.null(thin)&) {thin = round(0.9*niter/200)}
if(parameter=="Beta"){
chain=getBeta(fitted.model, burnin=burnin, thin=thin)
if(length(dim(chain)) == 3){
p=dim(chain)[2]
n.nets=dim(chain)[3]
}else{
p = 1
n.nets = dim(chain)[2] }
x=c(1:n.nets)
for(pie in 1:p){
int.quantiles=matrix(NA, n.nets,5)
for(kappa in 1:n.nets){
if(p > 1){
int.quantiles[kappa,]=quantile(chain[,pie,kappa],
probs=c(0.025, 0.25, 0.5, 0.75, 0.975))
}else{ int.quantiles[kappa,]=quantile(chain[,kappa],
probs=c(0.025, 0.25, 0.5, 0.75, 0.975)) }
}
plot(x, int.quantiles[,3], pch=".",
ylim=c(min(int.quantiles), max(int.quantiles)),
col="white", xlab="Network",
ylab=expression(beta[pie]),
main=paste("Beta", pie))
rect((x-0.05), int.quantiles[,1], (x+0.05),
int.quantiles[,5], col=colors[3], border=NA)
rect((x-0.15), int.quantiles[,2], (x+0.15),
int.quantiles[,4], col=colors[2], border=NA)
segments((x-0.15), int.quantiles[,3], (x+0.15),
int.quantiles[,3], lwd=2, col=colors[1])
}
}
if(parameter=="Intercept"){
chain=getIntercept(fitted.model, burnin=burnin, thin=thin)
p=1
n.nets=dim(chain)[2]
x=c(1:n.nets)
int.quantiles=matrix(NA, n.nets,5)
for(kappa in 1:n.nets){
int.quantiles[kappa,]=quantile(chain[,kappa],
probs=c(0.025, 0.25, 0.5, 0.75, 0.975))}
plot(x, int.quantiles[,3], pch=".",
ylim=c(min(int.quantiles), max(int.quantiles)),col="white",
xlab="Network", ylab=expression(beta[0]), main=expression(beta[0]))
rect((x-0.05), int.quantiles[,1], (x+0.05), int.quantiles[,5], col=colors[3], border=NA)
rect((x-0.15), int.quantiles[,2], (x+0.15), int.quantiles[,4], col=colors[2], border=NA)
segments((x-0.15), int.quantiles[,3], (x+0.15), int.quantiles[,3], lwd=2, col=colors[1])
}
if(parameter=="Theta"){
chain=getAlpha(fitted.model, burnin=burnin,thin=thin)
x=1
int.quantiles=matrix(NA, 1,5)
int.quantiles=quantile(chain, probs=c(0.025, 0.25, 0.5, 0.75, 0.975))
plot(x, int.quantiles[3], pch=".", ylim=c(min(int.quantiles),
max(int.quantiles)), col="white", xlab="Treatment Effect",
ylab=expression(theta), main=expression(theta))
rect((x-0.05), int.quantiles[1], (x+0.05), int.quantiles[5], col=colors[3], border=NA)
rect((x-0.15), int.quantiles[2], (x+0.15), int.quantiles[4], col=colors[2], border=NA)
segments((x-0.15), int.quantiles[3], (x+0.15), int.quantiles[3], lwd=2, col=colors[1])
}
}
###############
plotHLSM.fixed.fit=function(fitted.model, parameter, burnin =0, thin = 1){
colors=c("navy", "orange", "darkmagenta")
# if(is.null(burnin)){burnin = 0.1*fitted.model$call$niter}
# if(is.null(thin)) {thin = round(0.9*fitted.model$call$niter/200)}
if(parameter=="Beta"){
chain=getBeta(fitted.model,burnin=burnin,thin=thin)
if(is.null(dim(chain))==TRUE){p=1
}else {p=dim(chain)[2]}
x=c(1:p)
int.quantiles=matrix(NA, p,5)
if(p>1){
for(kappa in 1:p){
int.quantiles[kappa,]=quantile(chain[,kappa],
probs=c(0.025, 0.25, 0.5, 0.75, 0.975))}
}else{
kappa=1
int.quantiles[kappa,]=quantile(chain,
probs=c(0.025, 0.25, 0.5, 0.75, 0.975))}
plot(x, int.quantiles[,3], pch=".", lab=c(p, 5, 7),
xlim=c(0.75, p+0.25), ylim=c(min(int.quantiles), max(int.quantiles)),
col="white", xlab="Network", ylab=expression(beta), main=paste("Betas"))
rect((x-0.05), int.quantiles[,1], (x+0.05), int.quantiles[,5], col=colors[3], border=NA)
rect((x-0.15), int.quantiles[,2], (x+0.15), int.quantiles[,4], col=colors[2], border=NA)
segments((x-0.15), int.quantiles[,3], (x+0.15), int.quantiles[,3], lwd=2, col=colors[1])
}
if(parameter=="Intercept"){
chain=getIntercept(fitted.model,burnin=burnin, thin=thin)
x=1
int.quantiles=quantile(chain, probs=c(0.025, 0.25, 0.5, 0.75, 0.975))
plot(x, int.quantiles[3], pch=".", lab=c(1,5,7),
ylim=c(min(int.quantiles), max(int.quantiles)),
col="white", xlab="Network", ylab=expression(beta[0]),main=expression(beta[0]))
rect((x-0.05), int.quantiles[1], (x+0.05), int.quantiles[5], col=colors[3], border=NA)
rect((x-0.15), int.quantiles[2], (x+0.15), int.quantiles[4], col=colors[2], border=NA)
segments((x-0.15), int.quantiles[3], (x+0.15), int.quantiles[3], lwd=2, col=colors[1])
}
if(parameter=="Theta"){
chain=getAlpha(fitted.model,burnin=burnin,thin=thin)
x=1
int.quantiles=matrix(NA, 1,5)
int.quantiles=quantile(chain, probs=c(0.025, 0.25, 0.5, 0.75, 0.975))
plot(x, int.quantiles[3], pch=".", lab=c(1,5,7),
ylim=c(min(int.quantiles), max(int.quantiles)), col="white",
xlab="Treatment Effect", ylab=expression(theta), main=expression(theta))
rect((x-0.05), int.quantiles[1], (x+0.05), int.quantiles[5], col=colors[3], border=NA)
rect((x-0.15), int.quantiles[2], (x+0.15), int.quantiles[4], col=colors[2], border=NA)
segments((x-0.15), int.quantiles[3], (x+0.15), int.quantiles[3], lwd=2, col=colors[1])
}
}
######LATENT SPACE POSTIONS########
plotLS = function(LS,xx,fitted.model, node.name = FALSE,nodenames = NULL){
xcor = apply(LS[[xx]][,1,],1,mean)
ycor = apply(LS[[xx]][,2,],1,mean)
plot(xcor,ycor,main = paste('Network',xx),cex = 0.5,pch = 19,cex.lab = 1.5,col = 'red')
if(node.name == TRUE){
if(!is.null(nodenames)){
text(xcor,ycor,lab = nodenames) }
}
}
plotHLSM.LS = function(fitted.model,burnin=0,thin=1,...){
if(class(fitted.model)!='HLSM'){
stop('fitted.model must be of class HLSM')}
# niter = length(fitted.model$draws$Alpha)
# if(is.null(burnin)){burnin=0.1*niter}
# if(is.null(thin)){thin = round(0.9*niter/200)}
LS = getLS(fitted.model,burnin=burnin, thin=thin)
mat = matrix(1:length(LS),ceiling(length(LS)/5),5,byrow = TRUE)
dev.new(height = 10,width = 10)
layout(mat,widths = rep.int(1.5,ncol(mat)), heights = rep.int(1,nrow(mat)))
lapply(1:length(LS),function(x) plotLS(LS,x,fitted.model,...))
}