Nothing
# 3-d Tomography Plot for 2x3 Case
tomogRxC3d <- function(formula,data,total=NULL,lci=TRUE,estimates=FALSE,ci=FALSE,level=.95,seed=1234,color=hcl(h=30,c=100,l=60),transparency=.75,light=FALSE,rotate=TRUE){
##Run Through RxC Code Once (from Molly's original tomogRxC function)
#require(grDevices)
if (!requireNamespace("rgl", quietly = TRUE)) {
stop("Package rgl is needed for the tomogRxC3d function to work. Please install it.",
call. = FALSE)
}
noinfocount <- 0
form <- formula
dvname <- terms.formula(formula)[[2]]
covariate <- NA
#Make the bounds
rows <- c(all.names(form)[6:(length(all.names(form)))])
names <- rows
cols <- c(all.names(form)[3])
options(warn=-1)
bnds <- bounds(form, data=data,rows=rows, column =cols,threshold=0)
options(warn=0)
#Totals
dv <- data[, all.names(form)[3]]
#Assign other category
bndsoth <- bnds$bounds[[3]]
oth <- data[,all.names(form)[length(all.names(form))]]
#Assign x-axis category
bndsx <- bnds$bounds[[1]]
xcat <- data[,all.names(form)[6]]
#Assign y-axis category
bndsy <- bnds$bounds[[2]]
ycat <- data[,all.names(form)[7]]
#Minimums & Maximums
minx <- bndsx[,1]
miny <- bndsy[,1]
minoth <- bndsoth[,1]
maxx <- bndsx[,2]
maxy <- bndsy[,2]
maxoth <- bndsoth[,2]
#####
#Starting point when x is at minimum
##
#Holding x at its minimum, what are the bounds on y?
#When x is at its minimum, the new dv and total are:
newdv <- dv - (minx*xcat)
newtot <- oth + ycat
t <- newdv/newtot
y <- ycat/newtot
#The new bounds on the y category are:
lby <- cbind(miny, (t - maxoth*oth/newtot)/(y))
lby[,2] <- ifelse(y==0, 0, lby[,2])
lowy <- apply(lby,1,max)
hby <- cbind((t-minoth*oth/newtot)/y,maxy)
highy <- apply(hby,1,min)
#####
#Starting point when x is at maximum
##
#Holding x at its maximum, what are the bounds on y?
#The new bounds on x are:
newtot <- oth + xcat
newdv <- dv - (miny*ycat)
x <- xcat/newtot
t <- newdv/newtot
lbx <- cbind(minx, (t-maxoth*oth/newtot)/x)
lbx[,2] <- ifelse(x==0, 0, lbx[,2])
lowx <- apply(lbx,1,max)
hbx <- cbind((t-minoth*oth/newtot)/x,maxx)
highx <- apply(hbx,1,min)
#Graph starting points
#High starting points
hstr <- cbind(minx, highy)
#High ending points
hend <- cbind(highx, miny)
#Low starting points
lstr <- cbind(minx, lowy)
lend <- cbind(lowx, miny)
xl <- paste("Percent", names[1], dvname[2])
yl <- paste("Percent", names[2], dvname[2])
zl <- paste("Percent", names[3], dvname[2])
mn <- paste("3-d Tomography Plot in a 2x3 Table")
#plot(c(0,0), xlim=c(0,1), ylim=c(0,1), xaxs="i", yaxs="i",xlab=xl, ylab=yl, col="white", main=mn)
ok <- !is.na(hstr[,2]) & !is.na(hend[,1])
exp1 <- hstr[ok,2]>=maxy[ok]
exp2 <- hend[ok,1]>=maxx[ok]
exp3 <- lstr[ok,2]<=miny[ok]
exp4 <- lend[ok,1]<=minx[ok]
hstr <- hstr[ok,]
hend <- hend[ok,]
lstr <- lstr[ok,]
lend <- lend[ok,]
dv <- dv[ok]
ycat <- ycat[ok]
oth <- oth[ok]
minoth <- minoth[ok]
xcat <- xcat[ok]
maxy <- maxy[ok]
maxx <- maxx[ok]
##Create Matrices for coordinates
coords.xaxs<-matrix(data=NA,nrow=dim(hstr)[1],ncol=10)
coords.yaxs<-matrix(data=NA,nrow=dim(hstr)[1],ncol=10)
for(i in 1:dim(hstr)[1]){
if((exp1[i] + exp2[i] + exp3[i] + exp4[i])==0){
xaxs <- c(hstr[i,1], lstr[i,1],lend[i,1],hend[i,1])
yaxs <- c(hstr[i,2], lstr[i,2],lend[i,2], hend[i,2])
side1 <- c(xaxs,xaxs[1])
side2 <- c(yaxs,yaxs[1])
c1 <- sum(side1[1:(length(side1)-1)]*side2[2:length(side2)])
c2 <- sum(side1[2:(length(side1))]*side2[1:(length(side2)-1)])
area <- abs(c1-c2)/2
#polygon(xaxs,yaxs, col=rgb(0,0,0,alpha=0), border="black", lty=1, lwd=.25)
}
if((exp1[i]==1) & (exp2[i])==0){
cut <- (dv[i]-(oth[i])*minoth[i])/xcat[i] - maxy[i]*ycat[i]/xcat[i]
kink1x <- c(cut)
kink1y <- c(maxy[i])
}
if((exp2[i]==1) & (exp1[i])==0){
cut <- (dv[i]-(oth[i])*minoth[i])/ycat[i] - maxx[i]*xcat[i]/ycat[i]
kink1x <- c(maxx[i])
kink1y <- c(cut)
}
if((exp2[i]==1 & exp1[i]==1)){
cut <- (dv[i]-(oth[i])*minoth[i])/ycat[i] - maxx[i]*xcat[i]/ycat[i]
cut2 <- (dv[i]-(oth[i])*minoth[i])/xcat[i] - maxy[i]*ycat[i]/xcat[i]
kink1x <- c(maxx[i], cut2)
kink1y <- c(cut, maxy[i])
}
if((exp3[i]==1) & (exp4[i])==0){
cut <- (dv[i]-(oth[i])*maxoth[i])/xcat[i] - miny[i]*ycat[i]/xcat[i]
kink2x <- c(cut)
kink2y <- c(miny[i])
}
if((exp4[i]==1) & (exp3[i])==0){
cut <- (dv[i]-(oth[i])*maxoth[i])/ycat[i] - minx[i]*xcat[i]/ycat[i]
kink2x <- c(minx[i])
kink2y <- c(cut)
}
if((exp3[i]==1 & exp4[i]==1)){
cut <- (dv[i]-(oth[i])*maxoth[i])/ycat[i] - minx[i]*xcat[i]/ycat[i]
cut2 <- (dv[i]-(oth[i])*maxoth[i])/xcat[i] - miny[i]*ycat[i]/xcat[i]
kink2x <- c(minx[i], cut2)
kink2y <- c(cut, miny[i])
}
if((exp3[i] + exp4[i])==0 & (exp1[i] + exp2[i] + exp3[i] + exp4[i])!=0){
xaxs <- c(hstr[i,1], lstr[i,1],lend[i,1],hend[i,1], kink1x)
xaxs <- ifelse(is.nan(xaxs), 0, xaxs)
yaxs <- c(hstr[i,2], lstr[i,2],lend[i,2], hend[i,2], kink1y)
yaxs <- ifelse(is.nan(yaxs), 0, yaxs)
side1 <- c(xaxs,xaxs[1])
side2 <- c(yaxs,yaxs[1])
c1 <- sum(side1[1:(length(side1)-1)]*side2[2:length(side2)])
c2 <- sum(side1[2:(length(side1))]*side2[1:(length(side2)-1)])
area <- abs(c1-c2)/2
#polygon(xaxs,yaxs, col=rgb(0,0,0,alpha=0), border="black", lty=1, lwd=.25)
}
if((exp1[i] + exp2[i])==0 & (exp1[i] + exp2[i] + exp3[i] + exp4[i])!=0){
xaxs <- c(hstr[i,1], lstr[i,1],kink2x,lend[i,1],hend[i,1])
xaxs <- ifelse(is.nan(xaxs), 0, xaxs)
yaxs <- c(hstr[i,2], lstr[i,2],kink2y,lend[i,2], hend[i,2])
yaxs <- ifelse(is.nan(yaxs), 0, yaxs)
side1 <- c(xaxs,xaxs[1])
side2 <- c(yaxs,yaxs[1])
c1 <- sum(side1[1:(length(side1)-1)]*side2[2:length(side2)])
c2 <- sum(side1[2:(length(side1))]*side2[1:(length(side2)-1)])
area <- abs(c1-c2)/2
#polygon(xaxs,yaxs, col=rgb(0,0,0,alpha=0), border="black", lty=1, lwd=.25)
}
if((exp1[i] + exp2[i])!=0 & (exp3[i] + exp4[i])!=0){
xaxs <- c(hstr[i,1], lstr[i,1],kink2x,lend[i,1],hend[i,1], kink1x)
xaxs <- ifelse(is.nan(xaxs), 0, xaxs)
yaxs <- c(hstr[i,2], lstr[i,2],kink2y,lend[i,2], hend[i,2], kink1y)
yaxs <- ifelse(is.nan(yaxs), 0, yaxs)
side1 <- c(xaxs,xaxs[1])
side2 <- c(yaxs,yaxs[1])
c1 <- sum(side1[1:(length(side1)-1)]*side2[2:length(side2)])
c2 <- sum(side1[2:(length(side1))]*side2[1:(length(side2)-1)])
area <- abs(c1-c2)/2
#polygon(xaxs,yaxs, col=rgb(0,0,0,alpha=0), border="black", lty=1, lwd=.225)
}
for(j in 1:length(xaxs)){if(xaxs[j]<1e-5){xaxs[j]<-0}}
for(k in 1:length(yaxs)){if(yaxs[k]<1e-5){yaxs[k]<-0}}
for(l in 1:length(xaxs)){if(xaxs[l]>.99999){xaxs[l]<-1}}
for(m in 1:length(yaxs)){if(yaxs[m]>.99999){yaxs[m]<-1}}
axs<-rbind(xaxs,yaxs)
coords.axs<-unique(axs,MARGIN=2,incomparables=F)
coords.axs<-if(length(coords.axs[1,])!=10){cbind(coords.axs,matrix(data=NA,nrow=2,ncol=(10-length(coords.axs[1,]))))}else{coords.axs}
coords.xaxs[i,]<-coords.axs[1,]
coords.yaxs[i,]<-coords.axs[2,]
}
##Record Coordinate
w.coord<-coords.xaxs[,colSums(is.na(coords.xaxs))<nrow(coords.xaxs)]
b.coord<-coords.yaxs[,colSums(is.na(coords.yaxs))<nrow(coords.yaxs)]
#######################################################################
##THIS WILL BE A FN WHERE YOU SPECIFY NAME OF TURNOUT, TOTAL, CATEGORIES
##Combine Coordinates for 3d Graph
w<-matrix(NA,nrow=dim(hstr)[1],ncol=dim(w.coord)[2])
b<-matrix(NA,nrow=dim(hstr)[1],ncol=dim(w.coord)[2])
h<-matrix(NA,nrow=dim(hstr)[1],ncol=dim(w.coord)[2])
#Extract Data needed for accounting Identity
N <- data[,all.names(form)[3]] + data[,all.names(form)[4]]
t <- data[,all.names(form)[3]]/N
xb <- data[,all.names(form)[7]]/N
xw <- data[,all.names(form)[6]] /N
#Using the Accounting Identity from King (1997), calculate coordinates for remaining category
##T = B_b*X_b + B_w*X_w + B_h*X_h = B_b*X_b + B_w*X_w + B_h*(1-X_b-X_h)
coords<-matrix(NA, nrow=dim(w.coord)[2], ncol=3)
for(i in 1: dim(hstr)[1]){
coords<-cbind(w.coord[i,],b.coord[i,],rep(NA, nrow(coords)))
for(j in 1:nrow(coords)){
#Account for cases where xh=0
if(1-xb[i]-xw[i] < 1e-5){coords[j,3]<-0}else{
coords[j,3]<-(t[i] - coords[j,2]*xb[i] - coords[j,1]*xw[i])/(1-xb[i]-xw[i])}
##Make Sure that it records small numbers as 0
if(coords[j,3]<1e-5 & is.na(coords[j,3])==F){coords[j,3]<-0}
if(coords[j,3]>.99999 & is.na(coords[j,3])==F){coords[j,3]<-1}
##Account for cases where xb=0 and xw
if(xb[i]==0){coords[j,2]<-0}
if(xw[i]==0){coords[j,1]<-0}
}
w[i,]<-coords[,1]
b[i,]<-coords[,2]
h[i,]<-coords[,3]
}
##Run RxC Model
set.seed(seed)
dbuf <- ei(formula=formula, data=data)
##Calculate Center of Ellipsoid: mean of estimates
means<-apply(dbuf$draws$Beta, MARGIN=2, FUN=mean) #Mean of all simulations for each precinct and beta
means.w <- means[seq(1, length(means), 6)]
means.b <- means[seq(2, length(means), 6)]
means.h <- means[seq(3, length(means), 6)]
##Calculate Variance of Ellipsoid: average variance of simulations
ses.w <- sd(means.w)
ses.b <- sd(means.b)
ses.h <- sd(means.h)
##Calculate Covariance of Ellipsoid (Rho): average covariance of simulations
covs.bw<-cov(means.w,means.b)
covs.wh<-cov(means.w,means.h)
covs.bh<-cov(means.b,means.h)
#Create Rho Matrix
#Rho if covariance divided by multiplied variance
rho<-matrix(c(1,covs.bw/(ses.b*ses.w),covs.bh/(ses.b*ses.h),covs.bw/(ses.b*ses.w),1,covs.wh/(ses.h*ses.w)
,covs.bh/(ses.b*ses.h),covs.wh/(ses.h*ses.w),1),nrow=3)
######################################
##CALCULATE AREA OF PLANES for LoCoI
area <- rep(NA, nrow(w))
for(i in 1:nrow(w)){
#2-sided Line
if(length(unique(w[i,]))==1 | length(unique(b[i,]))==1 | length(unique(h[i,]))==1){
if(length(na.omit(w[i,]))>2){area[i]<-(sqrt(abs(h[i,1]-h[i,2])^2 + abs(b[i,1]-b[i,2])^2))*.01} #Multiply by epsilon to weight it appropriately small
if(length(na.omit(b[i,]))>2){area[i]<-(sqrt(abs(h[i,1]-h[i,2])^2 + abs(w[i,1]-w[i,2])^2))*.01}
if(length(na.omit(h[i,]))>2){area[i]<-(sqrt(abs(w[i,1]-w[i,2])^2 + abs(b[i,1]-b[i,2])^2))*.01}
}
# All other Planes
if(length(unique(w[i,]))!=1 & length(unique(b[i,]))!=1 & length(unique(h[i,]))!=1){
# Get 3d coordinates
p <- cbind(b[i,],w[i,],h[i,])
p <- na.omit(p)
# Get equation of plane
fit <- lm(h[i,] ~ b[i,] + w[i,])
coefs <- coef(fit)
a.eq <- coefs[2]
b.eq <- coefs[3]
c.eq <- -1
d.eq <- coefs["(Intercept)"]
# Get normal vector
n <- matrix(c(a.eq,b.eq,c.eq), nrow=1, ncol=3)
# Calculate Basis Vector on Plane
u1 <- matrix(c( (p[2,1]-p[1,1]), (p[2,2]-p[1,2]), (p[2,3]-p[1,3])), nrow=1, ncol=3)
# Get Second Basis Vector as cross-product of n and u1
CrossProduct3D <- function(x, y, i=1:3) {
# Project inputs into 3D, since the cross product only makes sense in 3D.
To3D <- function(x) head(c(x, rep(0, 3)), 3)
x <- To3D(x)
y <- To3D(y)
# Indices should be treated cyclically (i.e., index 4 is "really" index 1, and
# so on). Index3D() lets us do that using R's convention of 1-based (rather
# than 0-based) arrays.
Index3D <- function(i) (i - 1) %% 3 + 1
# The i'th component of the cross product is:
# (x[i + 1] * y[i + 2]) - (x[i + 2] * y[i + 1])
# as long as we treat the indices cyclically.
return (x[Index3D(i + 1)] * y[Index3D(i + 2)] -
x[Index3D(i + 2)] * y[Index3D(i + 1)])
}
u2 <- matrix(CrossProduct3D(n,u1), nrow=1, ncol=3)
# Normalize the basis vectors
e1 <- u1/sqrt(u1[1,1]^2 + u1[1,2]^2 + u1[1,3]^2)
e2 <- u2/sqrt(u2[1,1]^2 + u2[1,2]^2 + u2[1,3]^2)
# Get the projection matrix from basis vectors
pm <- rbind(e1, e2)
# Get the 2d coordinates from projection matrix
p.2d <- t(pm%*%t(p))
# Create matrix of points in counter-clockwise order (need to systematize)
ccorder <- rev(chull(p.2d)) #Convex hull gives clockwise order then reverse
parea <- p.2d[order(ccorder),]
# Add first row to bottom of matrix
parea <- rbind(parea, parea[1,])
# Calcluate area as 1/2 of determinant (need to systematize)
plus.area <- 0
for(j in 2:dim(parea)[1]){
x <- parea[j-1,1]*parea[j,2]
plus.area <- plus.area + x
}
minus.area <- 0
for(k in 2:dim(parea)[1]){
x <- parea[k-1,2]*parea[k,1]
minus.area <- minus.area + x
}
area[i] <- 1/2*(plus.area-minus.area)
}
}
scale<-(area-min(area))/(max(area)-min(area))*100
#################################
## Plot Based on Plot Type
##Tomography Plot - with/without point estimates
if(ci==F){
rgl::plot3d(b,w,h, type="s", col="black", size=0.3, xlim=c(0,1), ylim=c(0,1), zlim=c(0,1), xlab=yl, ylab=xl, zlab=zl)
##Add the plane for each set of points
for(i in 1:dim(hstr)[1]){
#Two point plane (line)
if(length(unique(w[i,]))==1 | length(unique(b[i,]))==1 | length(unique(h[i,]))==1){
if(lci==T){rgl::lines3d(x = b[i,], y = w[i,], z=h[i,], col=hcl(h=30,c=100,l=scale[i],alpha=1))}
else{rgl::lines3d(x = b[i,], y = w[i,], z=h[i,], col=color)}
}else{
#All Others
fit <- lm(h[i,] ~ b[i,] + w[i,])
coefs <- coef(fit)
a <- coefs[2]
c <- coefs[3]
d <- -1
e <- coefs["(Intercept)"]
if(lci==T){rgl::planes3d(a, c, d, e, alpha=transparency, col=hcl(h=30,c=100,l=scale[i]), lit=light)}
else{rgl::planes3d(a, c, d, e, alpha=transparency, col=color, lit=light)}
}
}
# With point estimates - mean of average simulations across precincts
if(estimates==T){rgl::spheres3d(means.b,means.w, means.h, col="black", radius=.03, add=T)}
}
##With Confidence Elipse
if(ci==T){
rgl::plot3d(rgl::ellipse3d(x=rho , scale= c(ses.b,ses.w,ses.h), centre= c(mean(means.b),mean(means.w), mean(means.h))
, alpha=.8, level=level) ,color=hcl(h=10,c=60,l=20), xlim=c(0,1), ylim=c(0,1), zlim=c(0,1), xlab=yl, ylab=xl, zlab=zl, lit=light)
rgl::points3d(mean(means.b),mean(means.w), mean(means.h), col="black", size=10, add=T) #Mean of average simulations across precincts
rgl::points3d(b,w,h, col="black", size=1, add=T)
##Add the plane for each set of points
for(i in 1:dim(hstr)[1]){
#Two point plane (line)
if(length(unique(w[i,]))==1 | length(unique(b[i,]))==1 | length(unique(h[i,]))==1){
if(lci==T){rgl::lines3d(x = b[i,], y = w[i,], z=h[i,], col=hcl(h=30,c=100,l=scale[i],alpha=1))}
else{rgl::lines3d(x = b[i,], y = w[i,], z=h[i,], col=hcl(h=30,c=100,l=60,alpha=1))}
}else{
#All Others
fit <- lm(h[i,] ~ b[i,] + w[i,])
coefs <- coef(fit)
a <- coefs[2]
c <- coefs[3]
d <- -1
e <- coefs["(Intercept)"]
if(lci==T){rgl::planes3d(a, c, d, e, alpha=.25, col=hcl(h=30,c=100,l=scale[i]), lit=light)}
else{rgl::planes3d(a, c, d, e, alpha=.25, col=hcl(h=30,c=100,l=60), lit=light)}
}
}
# With point estimates - mean of average simulations across precincts
if(estimates==T){rgl::spheres3d(means.b,means.w, means.h, col="black", radius=.03, add=T)}
}
if(rotate==T){rgl::play3d( rgl::spin3d(rpm=2.5), duration=20)}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.