###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=="Alpha"){ 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(alpha), main=expression(alpha)) 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=="Alpha"){ 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(alpha), main=expression(alpha)) 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,pdfname = NULL,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) if(!is.null(pdfname)){pdf(file=paste(pdfname,'.pdf',sep=''))}else(dev.new(height = 10,width = 10)) if(length(LS) > 5){ mat = matrix(1:length(LS),ceiling(length(LS)/5),5,byrow = TRUE) layout(mat,widths = rep.int(2.5,ncol(mat)), heights = rep.int(1,nrow(mat))) }else{ x1 = ceiling(length(LS)/2) par(mfrow = c(2,x1)) } lapply(1:length(LS),function(x) plotLS(LS,x,fitted.model,...)) if(!is.null(pdfname))dev.off() }