Nothing
R/Plot_jags.R
In IsotopeR: Stable Isotope Mixing Model
########output plots for methods paper###########
####Jake Ferguson 1/7/11
##updated 8/17/15
###########################
##Plots observations in 3d. Used with 3 isotopes
##
##########################
# for debugging
# > load('SampleOutput.Rdata')
# > source('Plot_jags.r'); RGL.plots(jags.out, X, sources, plot.mix=F, plot.ind.flag=FALSE, color.plots=TRUE)
#for diagnostics
# load('~/Dropbox/Work/workold/IsotopeR/R/SampleOutput.Rdata')
# source('Plot_jags.r')
# curves.plot(jags.1=jags.out, num.sources=num.sources, num.chains=mcmc.chains, color=T, individuals=N, xlab.vec=levels(as.factor(sources[,num.iso+1])), num.groups=num.groups)
# source('Plot_jags.r'); Tri.plots(jags.1=jags.out, X=X, sources=sources)
RGL.plots <- function(jags.1, X, sources, plot.mix=FALSE, plot.ind.flag=FALSE, color.plots=FALSE) {
requireNamespace("rgl", quietly = TRUE)
num.isos = 3 #this will always be true for these plots
sources.levs <- levels(sources[,num.isos+1])
num.sources <- nlevels(as.factor(sources.levs))
subsources <- sources[,num.isos+2]
groups = as.factor(X[,num.isos+1])
grays = grDevices::gray.colors(nlevels(groups), start = 0.3, end = 0.75, gamma = 2.2)
#define colors for groups
groupCols = vector('character', dim(X)[1])
if(nlevels(groups)==1) {groupCols = rep('dimgrey', dim(X)[1])} else {
for(i in 1:dim(X)[1]) {
groupCols[i] = grays[which(groups[i]==levels(groups))]
}
}
if(color.plots) {
##setting colors for the subsources
source.cols <- vector('numeric', dim(sources)[1])
source.list <- list()
source.subs <- list()
source.collist <- list()
h.vec <- seq(0,250, length.out=num.sources)
for(i in 1:num.sources) {
source.list[[i]] <- which(sources[,num.isos+1] == sources.levs[i])
source.subs[[i]] <- levels(as.factor(sources[source.list[[i]],num.isos+2]))
source.collist[[i]] <- colorspace::sequential_hcl(length(source.subs[[i]]), h = h.vec[i], c. = c(150, 10), l = c(30, 80), power = 1)
}
for(index in 1:dim(sources)[1]) {
source <- which(sources[index, num.isos+1] == sources.levs)
subsource <- which(as.factor(sources[index, num.isos+2]) == source.subs[[source]])
source.cols[index] <- source.collist[[source]][[subsource]]
}
} else {
source.cols <- vector('numeric')
source.cols= rep('black', dim(sources)[1])
}
if(!plot.mix) {
rgl::plot3d(x=sources[,1], y=sources[,2], z=sources[,3], xlab=names(sources)[1], ylab=names(sources)[2], zlab=names(sources)[3], col=source.cols, size=7, type="p")
rgl::plot3d(x=X[,1], y=X[,2], z=X[,3], size=10, col=groupCols, add=T)
rgl::decorate3d(xlab=NULL, ylab=NULL, zlab=NULL, main="Observations")
} else {
#plot the mixtures
jags.names <- dimnames(jags.1$summary$quantiles)
jags.table <- jags.1$summary$quantiles
num.ind <- length(unique(X[,5]))
med <- which(jags.names[[2]] == '50%')
low.ci <- which(jags.names[[2]] == '2.5%')
hi.ci <- which(jags.names[[2]] == '97.5%')
mix.index <- grep("mu.mix", jags.names[[1]])
cov.source.index <- grep("sd.source", jags.names[[1]])
mu.source.index <- grep("mu.source", jags.names[[1]])
mu.source <- jags.table[mu.source.index,med]
cov.source <- jags.table[cov.source.index, med]
mix.temp <- jags.table[mix.index,med]
mixmu.loCI <- jags.table[mix.index,low.ci]
mixmu.hiCI <- jags.table[mix.index,hi.ci]
mixmu.med <- cbind(mix.temp[1:num.ind], mix.temp[(num.ind+1):(2*num.ind)], mix.temp[(2*num.ind+1):(3*num.ind)])
mixmu.loCI <- cbind(mixmu.loCI[1:num.ind], mixmu.loCI[(num.ind+1):(2*num.ind)], mixmu.loCI[(2*num.ind+1):(3*num.ind)])
mixmu.hiCI <- cbind(mixmu.hiCI[1:num.ind], mixmu.hiCI[(num.ind+1):(2*num.ind)], mixmu.hiCI[(2*num.ind+1):(3*num.ind)])
colvec = vector('character', num.sources)
for(i in 1:num.sources) {
colvec[i] = source.collist[[i]][1]
}
rgl::plot3d(x=mu.source[1:num.sources], y=mu.source[1:num.sources+num.sources], z=mu.source[1:num.sources+2*num.sources], xlab=names(sources)[1], ylab=names(sources)[2], zlab=names(sources)[3], col=colvec, radius=0.25, type="s")
for(i in 1:num.sources) {
rgl::plot3d(x=mu.source[i] + c(-1,1)*cov.source[i], y=mu.source[num.sources+i], z=mu.source[2*num.sources+i] , type='l', add=T)
rgl::plot3d(x=mu.source[i], y=mu.source[num.sources+i] + c(-1,1)*cov.source[num.sources+i], z=mu.source[2*num.sources+i], type='l', add=T)
rgl::plot3d(x=mu.source[i], y=mu.source[num.sources+i], z=mu.source[2*num.sources+i] + c(-1,1)*cov.source[2*num.sources+i], type='l', add=T)
}
if(plot.ind.flag) {
rgl::plot3d(x=X[,1], y=X[,2], z=X[,3], size=10, col=groupCols, add=T)
} else {
rgl::plot3d(x=mixmu.med[,1], y=mixmu.med[,2], z=mixmu.med[,3], size=10, col=groupCols, add=T)
for(i in 1:dim(mixmu.med)[1]) {
rgl::plot3d(x=c(mixmu.loCI[i,1], mixmu.hiCI[i,1]), y=mixmu.med[i,2], z=mixmu.med[i,3] , type='l', add=T)
rgl::plot3d(x=mixmu.med[i,1], y=c(mixmu.loCI[i,2], mixmu.hiCI[i,2]), z=mixmu.med[i,3] , type='l', add=T)
rgl::plot3d(x=mixmu.med[i,1], y=mixmu.med[i,2], z=c(mixmu.loCI[i,3], mixmu.hiCI[i,3]) , type='l', add=T)
}
}
rgl::decorate3d(xlab=NULL, ylab=NULL, zlab=NULL, main="Estimated Mixing Space")
}
}
##################################
##Triangle plot in isotope space
##jags.1 is an jags object from an IsotopeR model run
##X is the observed isotope mixture values for each individual
##plot.mix is a 0,1 flag denoting whether to plot estimated mixture values
##plot.ind.flag is a 0,1 flag denoting whether to plot X.
##################################
Tri.plots <- function(jags.1, X, sources=NA, plot.mix=FALSE, plot.ind.flag=FALSE, me.flag=FALSE, color.plots=TRUE, xlab=NULL, ylab=NULL) {
requireNamespace("ellipse", quietly = TRUE)
num.isos <- 2 #this will always be true for these plots
num.ind <- dim(X)[1]
sources.levs <- levels(sources[,num.isos+1])
num.sources <- nlevels(as.factor(sources.levs))
subsources <- sources[,num.isos+2]
#groupC = 18+as.numeric(as.factor(X[,num.iso+1)]))
groups = as.factor(X[,num.isos+1])
grays = grDevices::gray.colors(nlevels(groups), start=0.3, end=0.75, gamma=2.2)
#define colors for groups
groupCols = vector('character', dim(X)[1])
if(nlevels(groups)==1) {groupCols = rep('dimgrey', dim(X)[1])} else {
for(i in 1:dim(X)[1]) {
groupCols[i] = grays[which(groups[i]==levels(groups))]
}
}
if(is.null(xlab) ) {
xlab=dimnames(X)[[2]][1]
}
if(is.null(ylab) ) {
ylab=dimnames(X)[[2]][2]
}
if(color.plots) {
##setting colors for the subsources
source.cols <- vector('numeric', dim(sources)[1])
source.list <- list()
source.subs <- list()
source.collist <- list()
h.vec <- seq(0,250, length.out=num.sources)
for(i in 1:num.sources) {
source.list[[i]] <- which(sources[,num.isos+1] == sources.levs[i])
source.subs[[i]] <- levels(as.factor(sources[source.list[[i]],num.isos+2]))
source.collist[[i]] <- colorspace::sequential_hcl(length(source.subs[[i]]), h = h.vec[i], c. = c(150, 10), l = c(30, 80), power = 1)
}
for(index in 1:dim(sources)[1]) {
source <- which(sources[index, num.isos+1] == sources.levs)
subsource <- which(as.factor(sources[index, num.isos+2]) == source.subs[[source]])
source.cols[index] <- source.collist[[source]][[subsource]]
}
} else {
source.cols <- vector('numeric')
for(levs in levels(sources[,3])) {source.cols[levs] <- "black" }
}
jags.names <- dimnames(jags.1$summary$quantiles)
jags.table <- jags.1$summary$quantiles
med <- which(jags.names[[2]] == '50%')
low.ci <- which(jags.names[[2]] == '2.5%')
hi.ci <- which(jags.names[[2]] == '97.5%')
mu.source.index <- grep("mu.source", jags.names[[1]])
mu.source <- jags.table[mu.source.index,med]
mu.lowCI <- jags.table[mu.source.index,low.ci]
mu.hiCI <- jags.table[mu.source.index,hi.ci]
cov.source.index <- grep("sd.source", jags.names[[1]])
cov.source <- jags.table[cov.source.index, med]
if(me.flag) {
sigmaz.index <- grep("sd.me", jags.names[[1]])
sigmaz.temp <- jags.table[sigmaz.index, med]
sigmaz.med <- sigmaz.temp
}
mix.index <- grep("mu.mix", jags.names[[1]])
mix.temp <- jags.table[mix.index,med]
mixmu.med <- cbind(mix.temp[1:num.ind], mix.temp[(num.ind+1):(2*num.ind)])
mixmu.loCI <- jags.table[mix.index,low.ci]
mixmu.loCI <- cbind(mixmu.loCI[1:num.ind], mixmu.loCI[(num.ind+1):(2*num.ind)])
mixmu.hiCI <- jags.table[mix.index,hi.ci]
mixmu.hiCI <- cbind(mixmu.hiCI[1:num.ind], mixmu.hiCI[(num.ind+1):(2*num.ind)])
cov.conc.index <- grep("sd.conc", jags.names[[1]])
cov.conc <- jags.table[cov.conc.index, med]
x.points <- mu.source[1:(length(mu.source)/2)]
x.sd <- cov.source[1:num.sources] #+ sigmaz.med[1]^2)
x.lowCI <- mu.lowCI[1:(length(mu.source)/2)]
x.hiCI <- mu.hiCI[1:(length(mu.source)/2)]
x.lowCI <- x.points - x.sd
x.hiCI <- x.points + x.sd
y.points <- mu.source[(length(mu.source)/2+1):length(mu.source)]
y.lowCI <- mu.lowCI[(length(mu.source)/2+1):length(mu.source)]
y.hiCI <- mu.hiCI[(length(mu.source)/2+1):length(mu.source)]
y.sd <- cov.source[(num.sources+1):(2*num.sources)]
y.lowCI <- y.points - 1*y.sd
y.hiCI <- y.points + 1*y.sd
data.cols="black"
D.source.index <- grep("mu.conc", jags.names[[1]])
D.source <- jags.table[D.source.index,med]
D.lowCI <- jags.table[D.source.index,low.ci]
D.hiCI <- jags.table[D.source.index,hi.ci]
if(length(D.source)==0) {
x.conc <- rep(1,3)
x.conc.sd <- rep(0,3)
y.conc <- rep(1,3)
y.conc.sd <- rep(0,3)
} else {
x.coord <- seq(1,length(D.source),by=2)
x.conc <- D.source[1:3]
x.conc.sd <- cov.source[1:3]
y.coord <- seq(2,length(D.source),by=2)
y.conc <- D.source[4:6]
y.conc.sd <- cov.source[4:6]
}
x.conc.lowCI <- x.conc - 1*x.conc.sd
x.conc.hiCI <- x.conc - 1*x.conc.sd
y.conc.lowCI <- y.conc - 1*y.conc.sd
y.conc.hiCI <- y.conc - 1*y.conc.sd
axis.raw <- seq(0,1,by=0.2)
base.matrix <- array(NA,c(length(axis.raw),length(axis.raw),3))
x.matrix <- base.matrix
xedge.matrix <- base.matrix
x.lowCI.matrix <- base.matrix
x.hiCI.matrix <- base.matrix
y.matrix <- base.matrix
yedge.matrix <- base.matrix
y.lowCI.matrix <- base.matrix
y.hiCI.matrix <- base.matrix
x.outer.matrix <- base.matrix
y.outer.matrix <- base.matrix
if(plot.mix & num.sources <=3) {
for(r1 in 1:length(axis.raw)) {
for(r2 in 1:(length(axis.raw))) {
temp <- 1- axis.raw[r1] - axis.raw[r2]
if(temp >= (0.00-1e-3) & temp <=1.001) {
base.matrix[r1,r2,1] <- axis.raw[r1]
base.matrix[r1,r2,2] <- axis.raw[r2]
base.matrix[r1,r2,3] <- 1- axis.raw[r1] - axis.raw[r2]
x.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*x.conc[1]/base.matrix[r1,r2,]%*%x.conc
x.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*x.conc[2]/base.matrix[r1,r2,]%*%x.conc
x.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*x.conc[3]/base.matrix[r1,r2,]%*%x.conc
y.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*y.conc[1]/base.matrix[r1,r2,]%*%y.conc
y.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*y.conc[2]/base.matrix[r1,r2,]%*%y.conc
y.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*y.conc[3]/base.matrix[r1,r2,]%*%y.conc
x.lowCI.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*x.conc.lowCI[1]/base.matrix[r1,r2,]%*%x.conc.lowCI
x.lowCI.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*x.conc.lowCI[2]/base.matrix[r1,r2,]%*%x.conc.lowCI
x.lowCI.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*x.conc.lowCI[3]/base.matrix[r1,r2,]%*%x.conc.lowCI
y.lowCI.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*y.conc.lowCI[1]/base.matrix[r1,r2,]%*%y.conc.lowCI
y.lowCI.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*y.conc.lowCI[2]/base.matrix[r1,r2,]%*%y.conc.lowCI
y.lowCI.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*y.conc.lowCI[3]/base.matrix[r1,r2,]%*%y.conc.lowCI
x.hiCI.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*x.conc.hiCI[1]/base.matrix[r1,r2,]%*%x.conc.hiCI
x.hiCI.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*x.conc.hiCI[2]/base.matrix[r1,r2,]%*%x.conc.hiCI
x.hiCI.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*x.conc.hiCI[3]/base.matrix[r1,r2,]%*%x.conc.hiCI
y.hiCI.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*y.conc.hiCI[1]/base.matrix[r1,r2,]%*%y.conc.hiCI
y.hiCI.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*y.conc.hiCI[2]/base.matrix[r1,r2,]%*%y.conc.hiCI
y.hiCI.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*y.conc.hiCI[3]/base.matrix[r1,r2,]%*%y.conc.hiCI
if(r1==1 | r2==1 | r1+r2 == (length(axis.raw)+1)) {
xedge.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*x.conc[1]/base.matrix[r1,r2,]%*%x.conc
xedge.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*x.conc[2]/base.matrix[r1,r2,]%*%x.conc
xedge.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*x.conc[3]/base.matrix[r1,r2,]%*%x.conc
yedge.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*y.conc[1]/base.matrix[r1,r2,]%*%y.conc
yedge.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*y.conc[2]/base.matrix[r1,r2,]%*%y.conc
yedge.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*y.conc[3]/base.matrix[r1,r2,]%*%y.conc
}
}
}
}#end r loop
x.basic <- apply(x.matrix, c(1,2), '%*%', x.points)
y.basic <- apply(y.matrix, c(1,2), '%*%', y.points)
x.edge <- apply(xedge.matrix, c(1,2), '%*%', x.points)
y.edge <- apply(yedge.matrix, c(1,2), '%*%', y.points)
x.edge.lowCI <- apply(x.lowCI.matrix, c(1,2), '%*%', x.points)
y.edge.lowCI <- apply(y.lowCI.matrix, c(1,2), '%*%', y.points)
x.edge.hiCI <- apply(x.hiCI.matrix, c(1,2), '%*%', x.points)
y.edge.hiCI <- apply(y.hiCI.matrix, c(1,2), '%*%', y.points)
x.outer.loCI <- apply(x.matrix, c(1,2), '%*%', x.lowCI)
x.outer.hiCI <- apply(x.matrix, c(1,2), '%*%', x.hiCI)
y.outer.loCI <- apply(y.matrix, c(1,2), '%*%', y.lowCI)
y.outer.hiCI <- apply(y.matrix, c(1,2), '%*%', y.hiCI)
}
if(!plot.ind.flag) {
plotrix::plotCI(x=x.points, y=y.points, liw=(x.sd), uiw=(x.sd), xlim=(range(x.points) + c(-1,1)*2.5*max(x.sd)), ylim=(range(y.points) + c(-1,1)*2.5*max(y.sd)), err="x",pch=19, xlab=xlab, ylab=ylab, col="white")
plotrix::plotCI(x=x.points, y=y.points, liw=(y.sd), uiw=(y.sd) ,err="y",add=TRUE,pch=19,col="white")
graphics::box(lwd=2)
base.matrix <- round(base.matrix,1)
if(num.sources <= 3 ) {
##draws the grey triangles
for(x1 in 1:(length(axis.raw))) {
for(y1 in 1:(length(axis.raw))) {
for(x2 in (x1-1):(x1+1)) {
for(y2 in (y1-1):(y1+1)) {
if(x2 > (length(axis.raw)) | x2 < 1.0 | y2 > (length(axis.raw)) | y2 < 1.0 ) {next()} else {
graphics::segments(x0=x.basic[x1,y1], y0=y.basic[x1,y1], x1=x.basic[x2,y2], y1=y.basic[x2,y2],col="grey",lwd=1,lty=1)
}
}
}
}#end for y1
}#end for x1
##draws the black edges & the CIs
edge.points <- which(!is.na(x.edge))
edge2 <- c(edge.points,edge.points)
back.vec <- (length(axis.raw)):1
for(i in 1:(length(axis.raw)-1)) {
##middle
graphics::segments(x0=x.basic[i,1], y0=y.basic[i,1], x1=x.basic[i+1,1], y1=y.basic[i+1,1],col="darkgrey",lwd=1.5)
graphics::segments(x0=x.basic[1,i], y0=y.basic[1,i], x1=x.basic[1,i+1], y1=y.basic[1,i+1],col="darkgrey",lwd=1.5)
graphics::segments(x0=x.basic[back.vec[i],i], y0=y.basic[back.vec[i],i], x1=x.basic[back.vec[i]-1,i+1], y1=y.basic[back.vec[i]-1,i+1],col="darkgrey",lwd=1.5)
}
}#end if num.sources
counter <- 1
for(levs in levels(sources[,3])) {
if(color.plots) { graphics::points(x.points[counter],y.points[counter],pch=14+counter,col=source.collist[[counter]],lwd=1) } else {
graphics::points(x.points[counter],y.points[counter],pch=14+counter,col="black",lwd=1) }
counter <- counter + 1
}
theta <- seq(0, 2 * pi, length=100)
rho.index <- grep("rho.mat", jags.names[[1]])
rho.vec <- jags.1$summary$statistics[rho.index,1]
rho.vec <- rho.vec[(num.sources+1):(2*num.sources)]
for(i in 1:num.sources) {
T.mat <- diag(2)
T.mat[1,2] <- rho.vec[i]
T.mat[2,1] <- rho.vec[i]
sd.vec <- c(cov.source[i], cov.source[i+3])
T.mat <- c(cov.source[i], cov.source[i], cov.source[i+3], cov.source[i+3])*T.mat*sd.vec
graphics::lines(ellipse::ellipse(T.mat, centre=c(x.points[i], y.points[i]), level=0.95), lty=2, lwd=2)
}
}#end plot.ind.flag
##plots predicted isotope values of individuals
if(plot.mix ==TRUE) {
plotrix::plotCI(x=mixmu.med[,1], y=mixmu.med[,2], liw=(mixmu.med[,1]-mixmu.loCI[,1]), uiw=(mixmu.hiCI[,1] - mixmu.med[,1]), sfrac=0, err="x", add=TRUE, col=groupCols, pch=19)
plotrix::plotCI(x=mixmu.med[,1], y=mixmu.med[,2], liw=(mixmu.med[,2]-mixmu.loCI[,2]), uiw=(mixmu.hiCI[,2] - mixmu.med[,2]), sfrac=0, err="y", , add=TRUE, col=groupCols, pch=19)
graphics::points(mixmu.med,lwd=1)
graphics::title("Estimated Mixing Space")
}
##Plots observed isotope values of individuals
if(plot.ind.flag) {
plotrix::plotCI(x=x.points, y=y.points, liw=(x.sd), uiw=(x.sd), sfrac=0, xlim=range(c(sources[,1], X[,1])), ylim=range(c(sources[,2], X[,2])), err="x",pch=19, xlab=xlab, ylab=ylab, col="white")
plotrix::plotCI(x=x.points, y=y.points, liw=(y.sd), uiw=(y.sd), sfrac=0, err="y",add=TRUE,pch=19,col="white")
box(lwd=2)
counter <- 0
for(levs in levels(sources[,3])) {
curr.lev <- which(sources[,3] == levs)
if(color.plots) {graphics::points(sources[curr.lev,1], sources[curr.lev,2], pch=15+counter, cex=1, col=source.cols[curr.lev])} else {
graphics::points(sources[curr.lev,1], sources[curr.lev,2], pch=15+counter, cex=1, col="black")}
counter <- counter+1
}
if(me.flag) {
plotrix::plotCI(x=X[,1], y=X[,2], uiw = 1*sqrt(sigmaz.med[1]), err="x", sfrac=0, add=TRUE, col=groupCols, pch=19)
plotrix::plotCI(x=X[,1], y=X[,2], uiw = 1*sqrt(sigmaz.med[2]), err="y", sfrac=0, add=TRUE, col=groupCols, pch=19)
}
graphics::points(X,pch=19,col=groupCols)
graphics::points(X,lwd=1)
graphics::title("Observations")
}
}#end triplots
##################################
##2 source plot in isotope space
##jags.1 is an jags object from an IsotopeR model run
##X is the observed isotope mixture values for each individual
##plot.mix is a 0,1 flag denoting whether to plot estimated mixture values
##plot.ind.flag is a 0,1 flag denoting whether to plot X
##################################
# source('R/Plot_jags.r')
Bi.plots <- function(jags.1, X, sources=NA, plot.mix=FALSE,plot.ind.flag=FALSE, me.flag=FALSE, color.plots=TRUE, xlab=NULL, ylab=NULL, xlim=NULL, ylim=NULL) {
N <- dim(X)[1]
num.isos <- 2#this should always be true for these plots
sources.levs <- levels(as.factor(sources[,num.isos+1]))
subsources <- sources[,num.isos+2]
source.cols <- vector('numeric', dim(sources)[1])
groups = as.factor(X[,num.isos+1])
grays = grDevices::gray.colors(nlevels(groups), start = 0.3, end = 0.75, gamma = 2.2)
#define colors for groups
groupCols = vector('character', dim(X)[1])
if(nlevels(groups)==1) {groupCols = rep('dimgrey', dim(X)[1])} else {
for(i in 1:dim(X)[1]) {
groupCols[i] = grays[which(groups[i]==levels(groups))]
}
}
if(color.plots) {
##setting colors for the subsources
source.list <- list()
source.subs <- list()
source.collist <- list()
source.list[[1]] <- which(as.factor(sources[,num.isos+1]) == sources.levs[1])
source.subs[[1]] <- levels(as.factor(sources[source.list[[1]],num.isos+2]))
source.collist[[1]] <- colorspace::sequential_hcl(length(source.subs[[1]]), h = 260, c. = c(150, 10), l = c(30, 80), power = 1)
source.list[[2]] <- which(as.factor(sources[,num.isos+1]) == sources.levs[2])
source.subs[[2]] <- levels(as.factor(sources[source.list[[2]],num.isos+2]))
source.collist[[2]] <- colorspace::sequential_hcl(length(source.subs[[2]]), h = 5, c. = c(200, 60), l = c(30, 90), power = 1)
for(index in 1:dim(sources)[1]) {
source <- which(as.factor(sources[index, num.isos+1]) == sources.levs)
subsource <- which(as.factor(sources[index, num.isos+2]) == source.subs[[source]])
source.cols[index] <- source.collist[[source]][[subsource]]
}
} else { source.cols <- "black" }
jags.names <- dimnames(jags.1$summary$quantiles)
jags.table <- jags.1$summary$quantiles
med <- which(jags.names[[2]] == '50%')
low.ci <- which(jags.names[[2]] == '2.5%')
hi.ci <- which(jags.names[[2]] == '97.5%')
mu.source.index <- grep("mu.source", jags.names[[1]])
mu.source <- jags.table[mu.source.index,med]
mu.lowCI <- jags.table[mu.source.index,low.ci]
mu.hiCI <- jags.table[mu.source.index,hi.ci]
cov.source.index <- grep("sd.source", jags.names[[1]])
cov.source <- jags.table[cov.source.index, med]
if(me.flag) {
sigmaz.index <- grep("sd.me", jags.names[[1]])
sigmaz.temp <- jags.table[sigmaz.index, med]
sigmaz.med <- sigmaz.temp#[c(1,4)]
}
mix.index <- grep("mu.mix", jags.names[[1]])
mix.temp <- jags.table[mix.index,med]
mixmu.med <- cbind(mix.temp[1:N], mix.temp[(N+1):(2*N)])
mixmu.loCI <- jags.table[mix.index,low.ci]
mixmu.loCI <- cbind(mixmu.loCI[1:N], mixmu.loCI[(N+1):(2*N)])
mixmu.hiCI <- jags.table[mix.index,hi.ci]
mixmu.hiCI <- cbind(mixmu.hiCI[1:N], mixmu.hiCI[(N+1):(2*N)])
cov.conc.index <- grep("sd.conc", jags.names[[1]])
cov.conc <- jags.table[cov.conc.index, med]
x.points <- mu.source[1:(length(mu.source)/2)]
x.sd <- cov.source[1:2] #+ sigmaz.med[1]^2)
x.lowCI <- x.points-1*x.sd
x.hiCI <- x.points+1*x.sd
y.points <- mu.source[(length(mu.source)/2+1):length(mu.source)]
y.lowCI <- mu.lowCI[(length(mu.source)/2+1):length(mu.source)]
y.hiCI <- mu.hiCI[(length(mu.source)/2+1):length(mu.source)]
y.sd <- cov.source[3:4]
y.lowCI <- y.points - 1*y.sd
y.hiCI <- y.points + 1*y.sd
data.cols="black"
D.source.index <- grep("mu.conc", jags.names[[1]])
D.source <- jags.table[D.source.index,med]
D.lowCI <- jags.table[D.source.index,low.ci]
D.hiCI <- jags.table[D.source.index,hi.ci]
if(length(D.source)==0) {
x.conc <- rep(1,2)
x.conc.sd <- rep(0,2)
y.conc <- rep(1,2)
y.conc.sd <- rep(0,2)
} else {
x.coord <- seq(1,length(D.source),by=2)
x.conc <- D.source[1:2]
x.conc.sd <- cov.source[1:2]
y.coord <- seq(2,length(D.source),by=2)
y.conc <- D.source[3:4]
y.conc.sd <- cov.source[3:4]
}
x.conc.lowCI <- x.conc-1*x.conc.sd
x.conc.hiCI <- x.conc+1*x.conc.sd
y.conc.lowCI <- y.conc-1*y.conc.sd
y.conc.hiCI <- y.conc+1*y.conc.sd
axis.raw <- seq(0,1,by=0.2)
base.matrix <- array(NA,c(length(axis.raw),2))
x.matrix <- base.matrix
xedge.matrix <- base.matrix
x.lowCI.matrix <- base.matrix
x.hiCI.matrix <- base.matrix
y.matrix <- base.matrix
yedge.matrix <- base.matrix
y.lowCI.matrix <- base.matrix
y.hiCI.matrix <- base.matrix
x.outer.matrix <- base.matrix
y.outer.matrix <- base.matrix
for(r1 in 1:length(axis.raw)) {
temp <- 1- axis.raw[r1] #- axis.raw[r2]
if(temp >= (0.00-1e-3) & temp <=1.001) {
base.matrix[r1,1] <- axis.raw[r1]
base.matrix[r1,2] <- 1-axis.raw[r1]#axis.raw[r2]
x.matrix[r1,1] <- base.matrix[r1,1]*x.conc[1]/base.matrix[r1,]%*%x.conc
x.matrix[r1,2] <- base.matrix[r1,2]*x.conc[2]/base.matrix[r1,]%*%x.conc
y.matrix[r1,1] <- base.matrix[r1,1]*y.conc[1]/base.matrix[r1,]%*%y.conc
y.matrix[r1,2] <- base.matrix[r1,2]*y.conc[2]/base.matrix[r1,]%*%y.conc
x.lowCI.matrix[r1,1] <- base.matrix[r1,1]*x.conc.lowCI[1]/base.matrix[r1,]%*%x.conc.lowCI
x.lowCI.matrix[r1,2] <- base.matrix[r1,2]*x.conc.lowCI[2]/base.matrix[r1,]%*%x.conc.lowCI
y.lowCI.matrix[r1,1] <- base.matrix[r1,1]*y.conc.lowCI[1]/base.matrix[r1,]%*%y.conc.lowCI
y.lowCI.matrix[r1,2] <- base.matrix[r1,2]*y.conc.lowCI[2]/base.matrix[r1,]%*%y.conc.lowCI
x.hiCI.matrix[r1,1] <- base.matrix[r1,1]*x.conc.hiCI[1]/base.matrix[r1,]%*%x.conc.hiCI
x.hiCI.matrix[r1,2] <- base.matrix[r1,2]*x.conc.hiCI[2]/base.matrix[r1,]%*%x.conc.hiCI
y.hiCI.matrix[r1,1] <- base.matrix[r1,1]*y.conc.hiCI[1]/base.matrix[r1,]%*%y.conc.hiCI
y.hiCI.matrix[r1,2] <- base.matrix[r1,2]*y.conc.hiCI[2]/base.matrix[r1,]%*%y.conc.hiCI
if(r1==1) {
xedge.matrix[r1,1] <- base.matrix[r1,1]*x.conc[1]/base.matrix[r1,]%*%x.conc
xedge.matrix[r1,2] <- base.matrix[r1,2]*x.conc[2]/base.matrix[r1,]%*%x.conc
yedge.matrix[r1,1] <- base.matrix[r1,1]*y.conc[1]/base.matrix[r1,]%*%y.conc
yedge.matrix[r1,2] <- base.matrix[r1,2]*y.conc[2]/base.matrix[r1,]%*%y.conc
}
}
}#end r loop
x.basic <- x.matrix%*%x.points
y.basic <- y.matrix%*%y.points
x.edge <- xedge.matrix%*%x.points
y.edge <- yedge.matrix%*%y.points
x.edge.lowCI <- x.lowCI.matrix%*%x.points
y.edge.lowCI <- y.lowCI.matrix%*%y.points
x.edge.hiCI <- x.hiCI.matrix%*%x.points
y.edge.hiCI <- y.hiCI.matrix%*%x.points
x.outer.loCI <- x.matrix%*%x.lowCI
x.outer.hiCI <- x.matrix%*%x.hiCI
y.outer.loCI <- y.matrix%*%y.lowCI
y.outer.hiCI <- y.matrix%*%y.hiCI
if( is.null(xlab) ) {
xlab=(dimnames(X)[[2]][1])
}
if( is.null(ylab) ) {
ylab=(dimnames(X)[[2]][2])
}
if(!plot.ind.flag) {
if(is.null(xlim)) { xlim <- (range(x.points) + c(-1,1)*3.0*x.sd) }
if(is.null(ylim)) { ylim <- range(y.points) + c(-1,1)*3.0*y.sd }
graphics::plot(x=x.points, y=y.points, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, pch=c(19,19), col=c("white","white"))
graphics::points(x.basic,y.basic,type='l',lwd=1,col="grey")
box(lwd=2)
if(plot.mix ==TRUE) {
plotrix::plotCI(x=mixmu.med[,1], y=mixmu.med[,2], liw=(mixmu.med[,1]-mixmu.loCI[,1]), uiw=(mixmu.hiCI[,1] - mixmu.med[,1]), sfrac=0, err="x", add=TRUE, col=c("black"),pch=19)
plotrix::plotCI(x=mixmu.med[,1], y=mixmu.med[,2], liw=(mixmu.med[,2]-mixmu.loCI[,2]), uiw=(mixmu.hiCI[,2] - mixmu.med[,2]), sfrac=0, err="y", add=TRUE, col=groupCols,pch=19)
graphics::points(mixmu.med,lwd=1)
if(color.plots) { graphics::points(x=x.points, y=y.points, pch=c(15,16),col=c(source.collist[[2]][1], source.collist[[1]][1])) } else { points(x=x.points, y=y.points, pch=c(15,16), col=c("black","black")) }
graphics::title('Estimated Mixing Space')
}
# rho.vec <- jags.1$BUGSoutput$mean$rho.source
rho.index <- grep("rho.mat", jags.names[[1]])
rho.vec <-jags.1$summary$statistics[rho.index,1]
rho.vec <- rho.vec[3:4]
for(i in 1:2) {
T.mat <- diag(2)
T.mat[1,2] <- rho.vec[i]
T.mat[2,1] <- rho.vec[i]
sd.vec <- c(cov.source[i], cov.source[i+2])
T.mat[1,2] <- rho.vec[i]^2
T.mat[2,1] <- rho.vec[i]^2
T.mat <- c(cov.source[i], cov.source[i], cov.source[i+2], cov.source[i+2])*T.mat*sd.vec
graphics::lines(ellipse::ellipse(T.mat, centre=c(x.points[i], y.points[i]), level=0.95), lty=2, lwd=2)
}
}#end plot
##Plots observed istope values of individuals
if(plot.ind.flag) {
if(is.null(xlim)) { xlim <- range(c(sources[,1],X[,1])) }
if(is.null(ylim)) { ylim <- range(c(sources[,2],X[,2])) }
plotrix::plotCI(x=x.points, y=y.points, liw=(x.sd), uiw=(x.sd), sfrac=0, xlim=xlim, ylim=ylim, err="x",pch=19, xlab=xlab, ylab=ylab, col="white")
plotrix::plotCI(x=x.points, y=y.points, liw=(y.sd), uiw=(y.sd), sfrac=0, err="y",add=TRUE,pch=19,col="white")
if(me.flag) {
plotrix::plotCI(x=X[,1], y=X[,2], uiw=1*sqrt(sigmaz.med[1]), sfrac=0, err="x", add=TRUE, col=groupCols, pch=19)
plotrix::plotCI(x=X[,1], y=X[,2], uiw=1*sqrt(sigmaz.med[2]), sfrac=0, err="y", add=TRUE, col=groupCols, pch=19)
} #else {
graphics::points(X, col=groupCols,pch=19)
graphics::points(X,lwd=1)
counter <- 0
for(levs in levels(as.factor(sources[,3]))) {
curr.lev <- which(as.factor(sources[,3]) == levs)
if(color.plots) {graphics::points(sources[curr.lev,1], sources[curr.lev,2], pch=15+counter, cex=1, col=source.cols[curr.lev])} else {
graphics::points(sources[curr.lev,1], sources[curr.lev,2], pch=15+counter, cex=1, col="black")}
counter <- counter+1
}
graphics::title("Observations")
}
}#end biplots
#######################################################
##plots smoothed histograms of individual distributions
#######################################################
curves.plot <- function(jags.1, num.sources, num.chains, color=FALSE, individuals, xlab.vec, num.groups=1) {
output <- jags.1$summary$quantiles
out.mcmc <- jags.1$mcmc
outnames.temp <- dimnames(jags.1$summary$quantiles)
out.names <- dimnames(output)[[1]]
out.perc <- which(outnames.temp[[2]] == "50%")
pIso.ind.array <- array(NA, c(2,individuals,3,3))
p.ind.array <- array(NA,c(individuals,3,3))
p1.popmed <- vector('numeric',num.sources)
mcmc.rbind <- NULL
for(i in 1:num.chains) {
mcmc.rbind <- rbind(mcmc.rbind, out.mcmc[[i]])
}
graphics::par(mfrow=c(num.sources,1))
source.colvec <- vector('numeric', num.sources)
h.vec <- seq(0,250, length.out=num.sources)
for(i in 1:num.sources) { source.colvec[i] <- colorspace::sequential_hcl(1, h = h.vec[i], c. = c(150, 10), l = c(30, 80), power = 1) }
pop.smooth.x <- matrix(NA, num.sources, 512)
pop.smooth.y <- matrix(NA, num.sources, 512)
ind.smooth.x <- list()
ind.smooth.y <- list()
group.smooth.x <- list()
group.smooth.y <- group.smooth.x
for(j in 1:num.sources) {
p.ind <- which(out.names == paste("p[",1,',',j,']',sep='') )
p.pop <- which(out.names == paste("p.pop[",j,']',sep='') )
if(num.groups > 1) { p.group <- which(out.names== paste("p.group[",j,']',sep='') ) }
p1.popmed[j] <- output[p.pop, out.perc]
temp <- stats::density(mcmc.rbind[,p.pop], kernel="epanechnikov", from=0, to=1)
pop.smooth.x[j,] <- temp$x
pop.smooth.y[j,] <- temp$y
ind.smooth.x[[j]] <- matrix(NA, individuals, 512)
ind.smooth.y[[j]] <- matrix(NA, individuals, 512)
if(length(p.ind) != 0) {
for(i in 1:individuals) {
p.ind <- which(out.names== paste("p[",i,',',j,']',sep='') )
temp <- density(mcmc.rbind[,p.ind], kernel="epanechnikov", from=0, to=1)
ind.smooth.x[[j]][i,] <- temp$x
ind.smooth.y[[j]][i,] <- temp$y
}
}
if(num.groups > 1) {
group.smooth.x[[j]] <- matrix(NA, num.groups, 512)
group.smooth.y[[j]] <- matrix(NA, num.groups, 512)
for(i in 1:num.groups) {
p.group <- which(out.names== paste("p.group[",i,',',j,']',sep='') )
temp.smooth <- density(mcmc.rbind[,p.group], kernel="epanechnikov", from=0, to=1)
group.smooth.x[[j]][i,] <- temp.smooth$x
group.smooth.y[[j]][i,] <- temp.smooth$y
}
}
}#end for
for(i in 1:num.sources) {
if(length(p.ind) == 0) {
if(length(group.smooth.y) != 0) {
ylims = c(0,max(pop.smooth.y[i,], group.smooth.y[[i]]))
} else {
ylims = c(0,max(pop.smooth.y[i,]))
}
graphics::plot(pop.smooth.x[i,], pop.smooth.y[i,], type='l', lwd=2, ylab="probability density", xlab=xlab.vec[j], main="", col=source.colvec[i])
}
else {
if(length(group.smooth.y) != 0) {
ylims = c(0,max(ind.smooth.y[[i]], pop.smooth.y[i,], group.smooth.y[[i]]))
} else {
ylims = c(0,max(ind.smooth.y[[i]], pop.smooth.y[i,]))
}
graphics::plot(pop.smooth.x[i,], pop.smooth.y[i,], type='l', lwd=2, xlim=c(0,1), ylim=ylims, ylab="probability density", xlab=xlab.vec[i], main="", col="white")
for(j in 1:individuals) {
if(color) { graphics::lines(ind.smooth.x[[i]][j,], ind.smooth.y[[i]][j,], type='l', lwd=2, col='grey') } else {
graphics::lines(ind.smooth.x[[i]][j,], ind.smooth.y[[i]][j,], type='l', lwd=2, col="grey") }
}
graphics::lines(pop.smooth.x[i,], pop.smooth.y[i,], type='l', lwd=2, ylab="probability density", xlab=xlab.vec[j], main="", col=source.colvec[i])
}
if(num.groups > 1) {
for(j in 1:num.groups) {
if(color) { graphics::lines(group.smooth.x[[i]][j,], group.smooth.y[[i]][j,], type='l', lwd=1, col='black', lty=2) } else {
graphics::lines(group.smooth.x[[i]][j,], group.smooth.y[[i]][j,], type='l', lwd=1, col="black", lty=2)}
}
}
}
}#end curves.plot
Try the IsotopeR package in your browser
Any scripts or data that you put into this service are public.
IsotopeR documentation built on May 2, 2019, 3:33 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.
########output plots for methods paper###########
####Jake Ferguson 1/7/11
##updated 8/17/15
###########################
##Plots observations in 3d. Used with 3 isotopes
##
##########################
# for debugging
# > load('SampleOutput.Rdata')
# > source('Plot_jags.r'); RGL.plots(jags.out, X, sources, plot.mix=F, plot.ind.flag=FALSE, color.plots=TRUE)
#for diagnostics
# load('~/Dropbox/Work/workold/IsotopeR/R/SampleOutput.Rdata')
# source('Plot_jags.r')
# curves.plot(jags.1=jags.out, num.sources=num.sources, num.chains=mcmc.chains, color=T, individuals=N, xlab.vec=levels(as.factor(sources[,num.iso+1])), num.groups=num.groups)
# source('Plot_jags.r'); Tri.plots(jags.1=jags.out, X=X, sources=sources)
RGL.plots <- function(jags.1, X, sources, plot.mix=FALSE, plot.ind.flag=FALSE, color.plots=FALSE) {
requireNamespace("rgl", quietly = TRUE)
num.isos = 3 #this will always be true for these plots
sources.levs <- levels(sources[,num.isos+1])
num.sources <- nlevels(as.factor(sources.levs))
subsources <- sources[,num.isos+2]
groups = as.factor(X[,num.isos+1])
grays = grDevices::gray.colors(nlevels(groups), start = 0.3, end = 0.75, gamma = 2.2)
#define colors for groups
groupCols = vector('character', dim(X)[1])
if(nlevels(groups)==1) {groupCols = rep('dimgrey', dim(X)[1])} else {
for(i in 1:dim(X)[1]) {
groupCols[i] = grays[which(groups[i]==levels(groups))]
}
}
if(color.plots) {
##setting colors for the subsources
source.cols <- vector('numeric', dim(sources)[1])
source.list <- list()
source.subs <- list()
source.collist <- list()
h.vec <- seq(0,250, length.out=num.sources)
for(i in 1:num.sources) {
source.list[[i]] <- which(sources[,num.isos+1] == sources.levs[i])
source.subs[[i]] <- levels(as.factor(sources[source.list[[i]],num.isos+2]))
source.collist[[i]] <- colorspace::sequential_hcl(length(source.subs[[i]]), h = h.vec[i], c. = c(150, 10), l = c(30, 80), power = 1)
}
for(index in 1:dim(sources)[1]) {
source <- which(sources[index, num.isos+1] == sources.levs)
subsource <- which(as.factor(sources[index, num.isos+2]) == source.subs[[source]])
source.cols[index] <- source.collist[[source]][[subsource]]
}
} else {
source.cols <- vector('numeric')
source.cols= rep('black', dim(sources)[1])
}
if(!plot.mix) {
rgl::plot3d(x=sources[,1], y=sources[,2], z=sources[,3], xlab=names(sources)[1], ylab=names(sources)[2], zlab=names(sources)[3], col=source.cols, size=7, type="p")
rgl::plot3d(x=X[,1], y=X[,2], z=X[,3], size=10, col=groupCols, add=T)
rgl::decorate3d(xlab=NULL, ylab=NULL, zlab=NULL, main="Observations")
} else {
#plot the mixtures
jags.names <- dimnames(jags.1$summary$quantiles)
jags.table <- jags.1$summary$quantiles
num.ind <- length(unique(X[,5]))
med <- which(jags.names[[2]] == '50%')
low.ci <- which(jags.names[[2]] == '2.5%')
hi.ci <- which(jags.names[[2]] == '97.5%')
mix.index <- grep("mu.mix", jags.names[[1]])
cov.source.index <- grep("sd.source", jags.names[[1]])
mu.source.index <- grep("mu.source", jags.names[[1]])
mu.source <- jags.table[mu.source.index,med]
cov.source <- jags.table[cov.source.index, med]
mix.temp <- jags.table[mix.index,med]
mixmu.loCI <- jags.table[mix.index,low.ci]
mixmu.hiCI <- jags.table[mix.index,hi.ci]
mixmu.med <- cbind(mix.temp[1:num.ind], mix.temp[(num.ind+1):(2*num.ind)], mix.temp[(2*num.ind+1):(3*num.ind)])
mixmu.loCI <- cbind(mixmu.loCI[1:num.ind], mixmu.loCI[(num.ind+1):(2*num.ind)], mixmu.loCI[(2*num.ind+1):(3*num.ind)])
mixmu.hiCI <- cbind(mixmu.hiCI[1:num.ind], mixmu.hiCI[(num.ind+1):(2*num.ind)], mixmu.hiCI[(2*num.ind+1):(3*num.ind)])
colvec = vector('character', num.sources)
for(i in 1:num.sources) {
colvec[i] = source.collist[[i]][1]
}
rgl::plot3d(x=mu.source[1:num.sources], y=mu.source[1:num.sources+num.sources], z=mu.source[1:num.sources+2*num.sources], xlab=names(sources)[1], ylab=names(sources)[2], zlab=names(sources)[3], col=colvec, radius=0.25, type="s")
for(i in 1:num.sources) {
rgl::plot3d(x=mu.source[i] + c(-1,1)*cov.source[i], y=mu.source[num.sources+i], z=mu.source[2*num.sources+i] , type='l', add=T)
rgl::plot3d(x=mu.source[i], y=mu.source[num.sources+i] + c(-1,1)*cov.source[num.sources+i], z=mu.source[2*num.sources+i], type='l', add=T)
rgl::plot3d(x=mu.source[i], y=mu.source[num.sources+i], z=mu.source[2*num.sources+i] + c(-1,1)*cov.source[2*num.sources+i], type='l', add=T)
}
if(plot.ind.flag) {
rgl::plot3d(x=X[,1], y=X[,2], z=X[,3], size=10, col=groupCols, add=T)
} else {
rgl::plot3d(x=mixmu.med[,1], y=mixmu.med[,2], z=mixmu.med[,3], size=10, col=groupCols, add=T)
for(i in 1:dim(mixmu.med)[1]) {
rgl::plot3d(x=c(mixmu.loCI[i,1], mixmu.hiCI[i,1]), y=mixmu.med[i,2], z=mixmu.med[i,3] , type='l', add=T)
rgl::plot3d(x=mixmu.med[i,1], y=c(mixmu.loCI[i,2], mixmu.hiCI[i,2]), z=mixmu.med[i,3] , type='l', add=T)
rgl::plot3d(x=mixmu.med[i,1], y=mixmu.med[i,2], z=c(mixmu.loCI[i,3], mixmu.hiCI[i,3]) , type='l', add=T)
}
}
rgl::decorate3d(xlab=NULL, ylab=NULL, zlab=NULL, main="Estimated Mixing Space")
}
}
##################################
##Triangle plot in isotope space
##jags.1 is an jags object from an IsotopeR model run
##X is the observed isotope mixture values for each individual
##plot.mix is a 0,1 flag denoting whether to plot estimated mixture values
##plot.ind.flag is a 0,1 flag denoting whether to plot X.
##################################
Tri.plots <- function(jags.1, X, sources=NA, plot.mix=FALSE, plot.ind.flag=FALSE, me.flag=FALSE, color.plots=TRUE, xlab=NULL, ylab=NULL) {
requireNamespace("ellipse", quietly = TRUE)
num.isos <- 2 #this will always be true for these plots
num.ind <- dim(X)[1]
sources.levs <- levels(sources[,num.isos+1])
num.sources <- nlevels(as.factor(sources.levs))
subsources <- sources[,num.isos+2]
#groupC = 18+as.numeric(as.factor(X[,num.iso+1)]))
groups = as.factor(X[,num.isos+1])
grays = grDevices::gray.colors(nlevels(groups), start=0.3, end=0.75, gamma=2.2)
#define colors for groups
groupCols = vector('character', dim(X)[1])
if(nlevels(groups)==1) {groupCols = rep('dimgrey', dim(X)[1])} else {
for(i in 1:dim(X)[1]) {
groupCols[i] = grays[which(groups[i]==levels(groups))]
}
}
if(is.null(xlab) ) {
xlab=dimnames(X)[[2]][1]
}
if(is.null(ylab) ) {
ylab=dimnames(X)[[2]][2]
}
if(color.plots) {
##setting colors for the subsources
source.cols <- vector('numeric', dim(sources)[1])
source.list <- list()
source.subs <- list()
source.collist <- list()
h.vec <- seq(0,250, length.out=num.sources)
for(i in 1:num.sources) {
source.list[[i]] <- which(sources[,num.isos+1] == sources.levs[i])
source.subs[[i]] <- levels(as.factor(sources[source.list[[i]],num.isos+2]))
source.collist[[i]] <- colorspace::sequential_hcl(length(source.subs[[i]]), h = h.vec[i], c. = c(150, 10), l = c(30, 80), power = 1)
}
for(index in 1:dim(sources)[1]) {
source <- which(sources[index, num.isos+1] == sources.levs)
subsource <- which(as.factor(sources[index, num.isos+2]) == source.subs[[source]])
source.cols[index] <- source.collist[[source]][[subsource]]
}
} else {
source.cols <- vector('numeric')
for(levs in levels(sources[,3])) {source.cols[levs] <- "black" }
}
jags.names <- dimnames(jags.1$summary$quantiles)
jags.table <- jags.1$summary$quantiles
med <- which(jags.names[[2]] == '50%')
low.ci <- which(jags.names[[2]] == '2.5%')
hi.ci <- which(jags.names[[2]] == '97.5%')
mu.source.index <- grep("mu.source", jags.names[[1]])
mu.source <- jags.table[mu.source.index,med]
mu.lowCI <- jags.table[mu.source.index,low.ci]
mu.hiCI <- jags.table[mu.source.index,hi.ci]
cov.source.index <- grep("sd.source", jags.names[[1]])
cov.source <- jags.table[cov.source.index, med]
if(me.flag) {
sigmaz.index <- grep("sd.me", jags.names[[1]])
sigmaz.temp <- jags.table[sigmaz.index, med]
sigmaz.med <- sigmaz.temp
}
mix.index <- grep("mu.mix", jags.names[[1]])
mix.temp <- jags.table[mix.index,med]
mixmu.med <- cbind(mix.temp[1:num.ind], mix.temp[(num.ind+1):(2*num.ind)])
mixmu.loCI <- jags.table[mix.index,low.ci]
mixmu.loCI <- cbind(mixmu.loCI[1:num.ind], mixmu.loCI[(num.ind+1):(2*num.ind)])
mixmu.hiCI <- jags.table[mix.index,hi.ci]
mixmu.hiCI <- cbind(mixmu.hiCI[1:num.ind], mixmu.hiCI[(num.ind+1):(2*num.ind)])
cov.conc.index <- grep("sd.conc", jags.names[[1]])
cov.conc <- jags.table[cov.conc.index, med]
x.points <- mu.source[1:(length(mu.source)/2)]
x.sd <- cov.source[1:num.sources] #+ sigmaz.med[1]^2)
x.lowCI <- mu.lowCI[1:(length(mu.source)/2)]
x.hiCI <- mu.hiCI[1:(length(mu.source)/2)]
x.lowCI <- x.points - x.sd
x.hiCI <- x.points + x.sd
y.points <- mu.source[(length(mu.source)/2+1):length(mu.source)]
y.lowCI <- mu.lowCI[(length(mu.source)/2+1):length(mu.source)]
y.hiCI <- mu.hiCI[(length(mu.source)/2+1):length(mu.source)]
y.sd <- cov.source[(num.sources+1):(2*num.sources)]
y.lowCI <- y.points - 1*y.sd
y.hiCI <- y.points + 1*y.sd
data.cols="black"
D.source.index <- grep("mu.conc", jags.names[[1]])
D.source <- jags.table[D.source.index,med]
D.lowCI <- jags.table[D.source.index,low.ci]
D.hiCI <- jags.table[D.source.index,hi.ci]
if(length(D.source)==0) {
x.conc <- rep(1,3)
x.conc.sd <- rep(0,3)
y.conc <- rep(1,3)
y.conc.sd <- rep(0,3)
} else {
x.coord <- seq(1,length(D.source),by=2)
x.conc <- D.source[1:3]
x.conc.sd <- cov.source[1:3]
y.coord <- seq(2,length(D.source),by=2)
y.conc <- D.source[4:6]
y.conc.sd <- cov.source[4:6]
}
x.conc.lowCI <- x.conc - 1*x.conc.sd
x.conc.hiCI <- x.conc - 1*x.conc.sd
y.conc.lowCI <- y.conc - 1*y.conc.sd
y.conc.hiCI <- y.conc - 1*y.conc.sd
axis.raw <- seq(0,1,by=0.2)
base.matrix <- array(NA,c(length(axis.raw),length(axis.raw),3))
x.matrix <- base.matrix
xedge.matrix <- base.matrix
x.lowCI.matrix <- base.matrix
x.hiCI.matrix <- base.matrix
y.matrix <- base.matrix
yedge.matrix <- base.matrix
y.lowCI.matrix <- base.matrix
y.hiCI.matrix <- base.matrix
x.outer.matrix <- base.matrix
y.outer.matrix <- base.matrix
if(plot.mix & num.sources <=3) {
for(r1 in 1:length(axis.raw)) {
for(r2 in 1:(length(axis.raw))) {
temp <- 1- axis.raw[r1] - axis.raw[r2]
if(temp >= (0.00-1e-3) & temp <=1.001) {
base.matrix[r1,r2,1] <- axis.raw[r1]
base.matrix[r1,r2,2] <- axis.raw[r2]
base.matrix[r1,r2,3] <- 1- axis.raw[r1] - axis.raw[r2]
x.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*x.conc[1]/base.matrix[r1,r2,]%*%x.conc
x.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*x.conc[2]/base.matrix[r1,r2,]%*%x.conc
x.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*x.conc[3]/base.matrix[r1,r2,]%*%x.conc
y.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*y.conc[1]/base.matrix[r1,r2,]%*%y.conc
y.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*y.conc[2]/base.matrix[r1,r2,]%*%y.conc
y.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*y.conc[3]/base.matrix[r1,r2,]%*%y.conc
x.lowCI.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*x.conc.lowCI[1]/base.matrix[r1,r2,]%*%x.conc.lowCI
x.lowCI.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*x.conc.lowCI[2]/base.matrix[r1,r2,]%*%x.conc.lowCI
x.lowCI.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*x.conc.lowCI[3]/base.matrix[r1,r2,]%*%x.conc.lowCI
y.lowCI.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*y.conc.lowCI[1]/base.matrix[r1,r2,]%*%y.conc.lowCI
y.lowCI.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*y.conc.lowCI[2]/base.matrix[r1,r2,]%*%y.conc.lowCI
y.lowCI.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*y.conc.lowCI[3]/base.matrix[r1,r2,]%*%y.conc.lowCI
x.hiCI.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*x.conc.hiCI[1]/base.matrix[r1,r2,]%*%x.conc.hiCI
x.hiCI.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*x.conc.hiCI[2]/base.matrix[r1,r2,]%*%x.conc.hiCI
x.hiCI.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*x.conc.hiCI[3]/base.matrix[r1,r2,]%*%x.conc.hiCI
y.hiCI.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*y.conc.hiCI[1]/base.matrix[r1,r2,]%*%y.conc.hiCI
y.hiCI.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*y.conc.hiCI[2]/base.matrix[r1,r2,]%*%y.conc.hiCI
y.hiCI.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*y.conc.hiCI[3]/base.matrix[r1,r2,]%*%y.conc.hiCI
if(r1==1 | r2==1 | r1+r2 == (length(axis.raw)+1)) {
xedge.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*x.conc[1]/base.matrix[r1,r2,]%*%x.conc
xedge.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*x.conc[2]/base.matrix[r1,r2,]%*%x.conc
xedge.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*x.conc[3]/base.matrix[r1,r2,]%*%x.conc
yedge.matrix[r1,r2,1] <- base.matrix[r1,r2,1]*y.conc[1]/base.matrix[r1,r2,]%*%y.conc
yedge.matrix[r1,r2,2] <- base.matrix[r1,r2,2]*y.conc[2]/base.matrix[r1,r2,]%*%y.conc
yedge.matrix[r1,r2,3] <- base.matrix[r1,r2,3]*y.conc[3]/base.matrix[r1,r2,]%*%y.conc
}
}
}
}#end r loop
x.basic <- apply(x.matrix, c(1,2), '%*%', x.points)
y.basic <- apply(y.matrix, c(1,2), '%*%', y.points)
x.edge <- apply(xedge.matrix, c(1,2), '%*%', x.points)
y.edge <- apply(yedge.matrix, c(1,2), '%*%', y.points)
x.edge.lowCI <- apply(x.lowCI.matrix, c(1,2), '%*%', x.points)
y.edge.lowCI <- apply(y.lowCI.matrix, c(1,2), '%*%', y.points)
x.edge.hiCI <- apply(x.hiCI.matrix, c(1,2), '%*%', x.points)
y.edge.hiCI <- apply(y.hiCI.matrix, c(1,2), '%*%', y.points)
x.outer.loCI <- apply(x.matrix, c(1,2), '%*%', x.lowCI)
x.outer.hiCI <- apply(x.matrix, c(1,2), '%*%', x.hiCI)
y.outer.loCI <- apply(y.matrix, c(1,2), '%*%', y.lowCI)
y.outer.hiCI <- apply(y.matrix, c(1,2), '%*%', y.hiCI)
}
if(!plot.ind.flag) {
plotrix::plotCI(x=x.points, y=y.points, liw=(x.sd), uiw=(x.sd), xlim=(range(x.points) + c(-1,1)*2.5*max(x.sd)), ylim=(range(y.points) + c(-1,1)*2.5*max(y.sd)), err="x",pch=19, xlab=xlab, ylab=ylab, col="white")
plotrix::plotCI(x=x.points, y=y.points, liw=(y.sd), uiw=(y.sd) ,err="y",add=TRUE,pch=19,col="white")
graphics::box(lwd=2)
base.matrix <- round(base.matrix,1)
if(num.sources <= 3 ) {
##draws the grey triangles
for(x1 in 1:(length(axis.raw))) {
for(y1 in 1:(length(axis.raw))) {
for(x2 in (x1-1):(x1+1)) {
for(y2 in (y1-1):(y1+1)) {
if(x2 > (length(axis.raw)) | x2 < 1.0 | y2 > (length(axis.raw)) | y2 < 1.0 ) {next()} else {
graphics::segments(x0=x.basic[x1,y1], y0=y.basic[x1,y1], x1=x.basic[x2,y2], y1=y.basic[x2,y2],col="grey",lwd=1,lty=1)
}
}
}
}#end for y1
}#end for x1
##draws the black edges & the CIs
edge.points <- which(!is.na(x.edge))
edge2 <- c(edge.points,edge.points)
back.vec <- (length(axis.raw)):1
for(i in 1:(length(axis.raw)-1)) {
##middle
graphics::segments(x0=x.basic[i,1], y0=y.basic[i,1], x1=x.basic[i+1,1], y1=y.basic[i+1,1],col="darkgrey",lwd=1.5)
graphics::segments(x0=x.basic[1,i], y0=y.basic[1,i], x1=x.basic[1,i+1], y1=y.basic[1,i+1],col="darkgrey",lwd=1.5)
graphics::segments(x0=x.basic[back.vec[i],i], y0=y.basic[back.vec[i],i], x1=x.basic[back.vec[i]-1,i+1], y1=y.basic[back.vec[i]-1,i+1],col="darkgrey",lwd=1.5)
}
}#end if num.sources
counter <- 1
for(levs in levels(sources[,3])) {
if(color.plots) { graphics::points(x.points[counter],y.points[counter],pch=14+counter,col=source.collist[[counter]],lwd=1) } else {
graphics::points(x.points[counter],y.points[counter],pch=14+counter,col="black",lwd=1) }
counter <- counter + 1
}
theta <- seq(0, 2 * pi, length=100)
rho.index <- grep("rho.mat", jags.names[[1]])
rho.vec <- jags.1$summary$statistics[rho.index,1]
rho.vec <- rho.vec[(num.sources+1):(2*num.sources)]
for(i in 1:num.sources) {
T.mat <- diag(2)
T.mat[1,2] <- rho.vec[i]
T.mat[2,1] <- rho.vec[i]
sd.vec <- c(cov.source[i], cov.source[i+3])
T.mat <- c(cov.source[i], cov.source[i], cov.source[i+3], cov.source[i+3])*T.mat*sd.vec
graphics::lines(ellipse::ellipse(T.mat, centre=c(x.points[i], y.points[i]), level=0.95), lty=2, lwd=2)
}
}#end plot.ind.flag
##plots predicted isotope values of individuals
if(plot.mix ==TRUE) {
plotrix::plotCI(x=mixmu.med[,1], y=mixmu.med[,2], liw=(mixmu.med[,1]-mixmu.loCI[,1]), uiw=(mixmu.hiCI[,1] - mixmu.med[,1]), sfrac=0, err="x", add=TRUE, col=groupCols, pch=19)
plotrix::plotCI(x=mixmu.med[,1], y=mixmu.med[,2], liw=(mixmu.med[,2]-mixmu.loCI[,2]), uiw=(mixmu.hiCI[,2] - mixmu.med[,2]), sfrac=0, err="y", , add=TRUE, col=groupCols, pch=19)
graphics::points(mixmu.med,lwd=1)
graphics::title("Estimated Mixing Space")
}
##Plots observed isotope values of individuals
if(plot.ind.flag) {
plotrix::plotCI(x=x.points, y=y.points, liw=(x.sd), uiw=(x.sd), sfrac=0, xlim=range(c(sources[,1], X[,1])), ylim=range(c(sources[,2], X[,2])), err="x",pch=19, xlab=xlab, ylab=ylab, col="white")
plotrix::plotCI(x=x.points, y=y.points, liw=(y.sd), uiw=(y.sd), sfrac=0, err="y",add=TRUE,pch=19,col="white")
box(lwd=2)
counter <- 0
for(levs in levels(sources[,3])) {
curr.lev <- which(sources[,3] == levs)
if(color.plots) {graphics::points(sources[curr.lev,1], sources[curr.lev,2], pch=15+counter, cex=1, col=source.cols[curr.lev])} else {
graphics::points(sources[curr.lev,1], sources[curr.lev,2], pch=15+counter, cex=1, col="black")}
counter <- counter+1
}
if(me.flag) {
plotrix::plotCI(x=X[,1], y=X[,2], uiw = 1*sqrt(sigmaz.med[1]), err="x", sfrac=0, add=TRUE, col=groupCols, pch=19)
plotrix::plotCI(x=X[,1], y=X[,2], uiw = 1*sqrt(sigmaz.med[2]), err="y", sfrac=0, add=TRUE, col=groupCols, pch=19)
}
graphics::points(X,pch=19,col=groupCols)
graphics::points(X,lwd=1)
graphics::title("Observations")
}
}#end triplots
##################################
##2 source plot in isotope space
##jags.1 is an jags object from an IsotopeR model run
##X is the observed isotope mixture values for each individual
##plot.mix is a 0,1 flag denoting whether to plot estimated mixture values
##plot.ind.flag is a 0,1 flag denoting whether to plot X
##################################
# source('R/Plot_jags.r')
Bi.plots <- function(jags.1, X, sources=NA, plot.mix=FALSE,plot.ind.flag=FALSE, me.flag=FALSE, color.plots=TRUE, xlab=NULL, ylab=NULL, xlim=NULL, ylim=NULL) {
N <- dim(X)[1]
num.isos <- 2#this should always be true for these plots
sources.levs <- levels(as.factor(sources[,num.isos+1]))
subsources <- sources[,num.isos+2]
source.cols <- vector('numeric', dim(sources)[1])
groups = as.factor(X[,num.isos+1])
grays = grDevices::gray.colors(nlevels(groups), start = 0.3, end = 0.75, gamma = 2.2)
#define colors for groups
groupCols = vector('character', dim(X)[1])
if(nlevels(groups)==1) {groupCols = rep('dimgrey', dim(X)[1])} else {
for(i in 1:dim(X)[1]) {
groupCols[i] = grays[which(groups[i]==levels(groups))]
}
}
if(color.plots) {
##setting colors for the subsources
source.list <- list()
source.subs <- list()
source.collist <- list()
source.list[[1]] <- which(as.factor(sources[,num.isos+1]) == sources.levs[1])
source.subs[[1]] <- levels(as.factor(sources[source.list[[1]],num.isos+2]))
source.collist[[1]] <- colorspace::sequential_hcl(length(source.subs[[1]]), h = 260, c. = c(150, 10), l = c(30, 80), power = 1)
source.list[[2]] <- which(as.factor(sources[,num.isos+1]) == sources.levs[2])
source.subs[[2]] <- levels(as.factor(sources[source.list[[2]],num.isos+2]))
source.collist[[2]] <- colorspace::sequential_hcl(length(source.subs[[2]]), h = 5, c. = c(200, 60), l = c(30, 90), power = 1)
for(index in 1:dim(sources)[1]) {
source <- which(as.factor(sources[index, num.isos+1]) == sources.levs)
subsource <- which(as.factor(sources[index, num.isos+2]) == source.subs[[source]])
source.cols[index] <- source.collist[[source]][[subsource]]
}
} else { source.cols <- "black" }
jags.names <- dimnames(jags.1$summary$quantiles)
jags.table <- jags.1$summary$quantiles
med <- which(jags.names[[2]] == '50%')
low.ci <- which(jags.names[[2]] == '2.5%')
hi.ci <- which(jags.names[[2]] == '97.5%')
mu.source.index <- grep("mu.source", jags.names[[1]])
mu.source <- jags.table[mu.source.index,med]
mu.lowCI <- jags.table[mu.source.index,low.ci]
mu.hiCI <- jags.table[mu.source.index,hi.ci]
cov.source.index <- grep("sd.source", jags.names[[1]])
cov.source <- jags.table[cov.source.index, med]
if(me.flag) {
sigmaz.index <- grep("sd.me", jags.names[[1]])
sigmaz.temp <- jags.table[sigmaz.index, med]
sigmaz.med <- sigmaz.temp#[c(1,4)]
}
mix.index <- grep("mu.mix", jags.names[[1]])
mix.temp <- jags.table[mix.index,med]
mixmu.med <- cbind(mix.temp[1:N], mix.temp[(N+1):(2*N)])
mixmu.loCI <- jags.table[mix.index,low.ci]
mixmu.loCI <- cbind(mixmu.loCI[1:N], mixmu.loCI[(N+1):(2*N)])
mixmu.hiCI <- jags.table[mix.index,hi.ci]
mixmu.hiCI <- cbind(mixmu.hiCI[1:N], mixmu.hiCI[(N+1):(2*N)])
cov.conc.index <- grep("sd.conc", jags.names[[1]])
cov.conc <- jags.table[cov.conc.index, med]
x.points <- mu.source[1:(length(mu.source)/2)]
x.sd <- cov.source[1:2] #+ sigmaz.med[1]^2)
x.lowCI <- x.points-1*x.sd
x.hiCI <- x.points+1*x.sd
y.points <- mu.source[(length(mu.source)/2+1):length(mu.source)]
y.lowCI <- mu.lowCI[(length(mu.source)/2+1):length(mu.source)]
y.hiCI <- mu.hiCI[(length(mu.source)/2+1):length(mu.source)]
y.sd <- cov.source[3:4]
y.lowCI <- y.points - 1*y.sd
y.hiCI <- y.points + 1*y.sd
data.cols="black"
D.source.index <- grep("mu.conc", jags.names[[1]])
D.source <- jags.table[D.source.index,med]
D.lowCI <- jags.table[D.source.index,low.ci]
D.hiCI <- jags.table[D.source.index,hi.ci]
if(length(D.source)==0) {
x.conc <- rep(1,2)
x.conc.sd <- rep(0,2)
y.conc <- rep(1,2)
y.conc.sd <- rep(0,2)
} else {
x.coord <- seq(1,length(D.source),by=2)
x.conc <- D.source[1:2]
x.conc.sd <- cov.source[1:2]
y.coord <- seq(2,length(D.source),by=2)
y.conc <- D.source[3:4]
y.conc.sd <- cov.source[3:4]
}
x.conc.lowCI <- x.conc-1*x.conc.sd
x.conc.hiCI <- x.conc+1*x.conc.sd
y.conc.lowCI <- y.conc-1*y.conc.sd
y.conc.hiCI <- y.conc+1*y.conc.sd
axis.raw <- seq(0,1,by=0.2)
base.matrix <- array(NA,c(length(axis.raw),2))
x.matrix <- base.matrix
xedge.matrix <- base.matrix
x.lowCI.matrix <- base.matrix
x.hiCI.matrix <- base.matrix
y.matrix <- base.matrix
yedge.matrix <- base.matrix
y.lowCI.matrix <- base.matrix
y.hiCI.matrix <- base.matrix
x.outer.matrix <- base.matrix
y.outer.matrix <- base.matrix
for(r1 in 1:length(axis.raw)) {
temp <- 1- axis.raw[r1] #- axis.raw[r2]
if(temp >= (0.00-1e-3) & temp <=1.001) {
base.matrix[r1,1] <- axis.raw[r1]
base.matrix[r1,2] <- 1-axis.raw[r1]#axis.raw[r2]
x.matrix[r1,1] <- base.matrix[r1,1]*x.conc[1]/base.matrix[r1,]%*%x.conc
x.matrix[r1,2] <- base.matrix[r1,2]*x.conc[2]/base.matrix[r1,]%*%x.conc
y.matrix[r1,1] <- base.matrix[r1,1]*y.conc[1]/base.matrix[r1,]%*%y.conc
y.matrix[r1,2] <- base.matrix[r1,2]*y.conc[2]/base.matrix[r1,]%*%y.conc
x.lowCI.matrix[r1,1] <- base.matrix[r1,1]*x.conc.lowCI[1]/base.matrix[r1,]%*%x.conc.lowCI
x.lowCI.matrix[r1,2] <- base.matrix[r1,2]*x.conc.lowCI[2]/base.matrix[r1,]%*%x.conc.lowCI
y.lowCI.matrix[r1,1] <- base.matrix[r1,1]*y.conc.lowCI[1]/base.matrix[r1,]%*%y.conc.lowCI
y.lowCI.matrix[r1,2] <- base.matrix[r1,2]*y.conc.lowCI[2]/base.matrix[r1,]%*%y.conc.lowCI
x.hiCI.matrix[r1,1] <- base.matrix[r1,1]*x.conc.hiCI[1]/base.matrix[r1,]%*%x.conc.hiCI
x.hiCI.matrix[r1,2] <- base.matrix[r1,2]*x.conc.hiCI[2]/base.matrix[r1,]%*%x.conc.hiCI
y.hiCI.matrix[r1,1] <- base.matrix[r1,1]*y.conc.hiCI[1]/base.matrix[r1,]%*%y.conc.hiCI
y.hiCI.matrix[r1,2] <- base.matrix[r1,2]*y.conc.hiCI[2]/base.matrix[r1,]%*%y.conc.hiCI
if(r1==1) {
xedge.matrix[r1,1] <- base.matrix[r1,1]*x.conc[1]/base.matrix[r1,]%*%x.conc
xedge.matrix[r1,2] <- base.matrix[r1,2]*x.conc[2]/base.matrix[r1,]%*%x.conc
yedge.matrix[r1,1] <- base.matrix[r1,1]*y.conc[1]/base.matrix[r1,]%*%y.conc
yedge.matrix[r1,2] <- base.matrix[r1,2]*y.conc[2]/base.matrix[r1,]%*%y.conc
}
}
}#end r loop
x.basic <- x.matrix%*%x.points
y.basic <- y.matrix%*%y.points
x.edge <- xedge.matrix%*%x.points
y.edge <- yedge.matrix%*%y.points
x.edge.lowCI <- x.lowCI.matrix%*%x.points
y.edge.lowCI <- y.lowCI.matrix%*%y.points
x.edge.hiCI <- x.hiCI.matrix%*%x.points
y.edge.hiCI <- y.hiCI.matrix%*%x.points
x.outer.loCI <- x.matrix%*%x.lowCI
x.outer.hiCI <- x.matrix%*%x.hiCI
y.outer.loCI <- y.matrix%*%y.lowCI
y.outer.hiCI <- y.matrix%*%y.hiCI
if( is.null(xlab) ) {
xlab=(dimnames(X)[[2]][1])
}
if( is.null(ylab) ) {
ylab=(dimnames(X)[[2]][2])
}
if(!plot.ind.flag) {
if(is.null(xlim)) { xlim <- (range(x.points) + c(-1,1)*3.0*x.sd) }
if(is.null(ylim)) { ylim <- range(y.points) + c(-1,1)*3.0*y.sd }
graphics::plot(x=x.points, y=y.points, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, pch=c(19,19), col=c("white","white"))
graphics::points(x.basic,y.basic,type='l',lwd=1,col="grey")
box(lwd=2)
if(plot.mix ==TRUE) {
plotrix::plotCI(x=mixmu.med[,1], y=mixmu.med[,2], liw=(mixmu.med[,1]-mixmu.loCI[,1]), uiw=(mixmu.hiCI[,1] - mixmu.med[,1]), sfrac=0, err="x", add=TRUE, col=c("black"),pch=19)
plotrix::plotCI(x=mixmu.med[,1], y=mixmu.med[,2], liw=(mixmu.med[,2]-mixmu.loCI[,2]), uiw=(mixmu.hiCI[,2] - mixmu.med[,2]), sfrac=0, err="y", add=TRUE, col=groupCols,pch=19)
graphics::points(mixmu.med,lwd=1)
if(color.plots) { graphics::points(x=x.points, y=y.points, pch=c(15,16),col=c(source.collist[[2]][1], source.collist[[1]][1])) } else { points(x=x.points, y=y.points, pch=c(15,16), col=c("black","black")) }
graphics::title('Estimated Mixing Space')
}
# rho.vec <- jags.1$BUGSoutput$mean$rho.source
rho.index <- grep("rho.mat", jags.names[[1]])
rho.vec <-jags.1$summary$statistics[rho.index,1]
rho.vec <- rho.vec[3:4]
for(i in 1:2) {
T.mat <- diag(2)
T.mat[1,2] <- rho.vec[i]
T.mat[2,1] <- rho.vec[i]
sd.vec <- c(cov.source[i], cov.source[i+2])
T.mat[1,2] <- rho.vec[i]^2
T.mat[2,1] <- rho.vec[i]^2
T.mat <- c(cov.source[i], cov.source[i], cov.source[i+2], cov.source[i+2])*T.mat*sd.vec
graphics::lines(ellipse::ellipse(T.mat, centre=c(x.points[i], y.points[i]), level=0.95), lty=2, lwd=2)
}
}#end plot
##Plots observed istope values of individuals
if(plot.ind.flag) {
if(is.null(xlim)) { xlim <- range(c(sources[,1],X[,1])) }
if(is.null(ylim)) { ylim <- range(c(sources[,2],X[,2])) }
plotrix::plotCI(x=x.points, y=y.points, liw=(x.sd), uiw=(x.sd), sfrac=0, xlim=xlim, ylim=ylim, err="x",pch=19, xlab=xlab, ylab=ylab, col="white")
plotrix::plotCI(x=x.points, y=y.points, liw=(y.sd), uiw=(y.sd), sfrac=0, err="y",add=TRUE,pch=19,col="white")
if(me.flag) {
plotrix::plotCI(x=X[,1], y=X[,2], uiw=1*sqrt(sigmaz.med[1]), sfrac=0, err="x", add=TRUE, col=groupCols, pch=19)
plotrix::plotCI(x=X[,1], y=X[,2], uiw=1*sqrt(sigmaz.med[2]), sfrac=0, err="y", add=TRUE, col=groupCols, pch=19)
} #else {
graphics::points(X, col=groupCols,pch=19)
graphics::points(X,lwd=1)
counter <- 0
for(levs in levels(as.factor(sources[,3]))) {
curr.lev <- which(as.factor(sources[,3]) == levs)
if(color.plots) {graphics::points(sources[curr.lev,1], sources[curr.lev,2], pch=15+counter, cex=1, col=source.cols[curr.lev])} else {
graphics::points(sources[curr.lev,1], sources[curr.lev,2], pch=15+counter, cex=1, col="black")}
counter <- counter+1
}
graphics::title("Observations")
}
}#end biplots
#######################################################
##plots smoothed histograms of individual distributions
#######################################################
curves.plot <- function(jags.1, num.sources, num.chains, color=FALSE, individuals, xlab.vec, num.groups=1) {
output <- jags.1$summary$quantiles
out.mcmc <- jags.1$mcmc
outnames.temp <- dimnames(jags.1$summary$quantiles)
out.names <- dimnames(output)[[1]]
out.perc <- which(outnames.temp[[2]] == "50%")
pIso.ind.array <- array(NA, c(2,individuals,3,3))
p.ind.array <- array(NA,c(individuals,3,3))
p1.popmed <- vector('numeric',num.sources)
mcmc.rbind <- NULL
for(i in 1:num.chains) {
mcmc.rbind <- rbind(mcmc.rbind, out.mcmc[[i]])
}
graphics::par(mfrow=c(num.sources,1))
source.colvec <- vector('numeric', num.sources)
h.vec <- seq(0,250, length.out=num.sources)
for(i in 1:num.sources) { source.colvec[i] <- colorspace::sequential_hcl(1, h = h.vec[i], c. = c(150, 10), l = c(30, 80), power = 1) }
pop.smooth.x <- matrix(NA, num.sources, 512)
pop.smooth.y <- matrix(NA, num.sources, 512)
ind.smooth.x <- list()
ind.smooth.y <- list()
group.smooth.x <- list()
group.smooth.y <- group.smooth.x
for(j in 1:num.sources) {
p.ind <- which(out.names == paste("p[",1,',',j,']',sep='') )
p.pop <- which(out.names == paste("p.pop[",j,']',sep='') )
if(num.groups > 1) { p.group <- which(out.names== paste("p.group[",j,']',sep='') ) }
p1.popmed[j] <- output[p.pop, out.perc]
temp <- stats::density(mcmc.rbind[,p.pop], kernel="epanechnikov", from=0, to=1)
pop.smooth.x[j,] <- temp$x
pop.smooth.y[j,] <- temp$y
ind.smooth.x[[j]] <- matrix(NA, individuals, 512)
ind.smooth.y[[j]] <- matrix(NA, individuals, 512)
if(length(p.ind) != 0) {
for(i in 1:individuals) {
p.ind <- which(out.names== paste("p[",i,',',j,']',sep='') )
temp <- density(mcmc.rbind[,p.ind], kernel="epanechnikov", from=0, to=1)
ind.smooth.x[[j]][i,] <- temp$x
ind.smooth.y[[j]][i,] <- temp$y
}
}
if(num.groups > 1) {
group.smooth.x[[j]] <- matrix(NA, num.groups, 512)
group.smooth.y[[j]] <- matrix(NA, num.groups, 512)
for(i in 1:num.groups) {
p.group <- which(out.names== paste("p.group[",i,',',j,']',sep='') )
temp.smooth <- density(mcmc.rbind[,p.group], kernel="epanechnikov", from=0, to=1)
group.smooth.x[[j]][i,] <- temp.smooth$x
group.smooth.y[[j]][i,] <- temp.smooth$y
}
}
}#end for
for(i in 1:num.sources) {
if(length(p.ind) == 0) {
if(length(group.smooth.y) != 0) {
ylims = c(0,max(pop.smooth.y[i,], group.smooth.y[[i]]))
} else {
ylims = c(0,max(pop.smooth.y[i,]))
}
graphics::plot(pop.smooth.x[i,], pop.smooth.y[i,], type='l', lwd=2, ylab="probability density", xlab=xlab.vec[j], main="", col=source.colvec[i])
}
else {
if(length(group.smooth.y) != 0) {
ylims = c(0,max(ind.smooth.y[[i]], pop.smooth.y[i,], group.smooth.y[[i]]))
} else {
ylims = c(0,max(ind.smooth.y[[i]], pop.smooth.y[i,]))
}
graphics::plot(pop.smooth.x[i,], pop.smooth.y[i,], type='l', lwd=2, xlim=c(0,1), ylim=ylims, ylab="probability density", xlab=xlab.vec[i], main="", col="white")
for(j in 1:individuals) {
if(color) { graphics::lines(ind.smooth.x[[i]][j,], ind.smooth.y[[i]][j,], type='l', lwd=2, col='grey') } else {
graphics::lines(ind.smooth.x[[i]][j,], ind.smooth.y[[i]][j,], type='l', lwd=2, col="grey") }
}
graphics::lines(pop.smooth.x[i,], pop.smooth.y[i,], type='l', lwd=2, ylab="probability density", xlab=xlab.vec[j], main="", col=source.colvec[i])
}
if(num.groups > 1) {
for(j in 1:num.groups) {
if(color) { graphics::lines(group.smooth.x[[i]][j,], group.smooth.y[[i]][j,], type='l', lwd=1, col='black', lty=2) } else {
graphics::lines(group.smooth.x[[i]][j,], group.smooth.y[[i]][j,], type='l', lwd=1, col="black", lty=2)}
}
}
}
}#end curves.plot
Try the IsotopeR 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.