################################################################################### # # Contingency table and 3D Histogram of two variables # # Last review: 17/06/2018 # email: ja@fct.unl.pt # ################################################################################### ################################################################################### varlistindex <- c(5,6) nclasses <- c(10,12) titleUNITS <- c("ppm","ppm") spaceBARS <- 0.5 angles <- c(45,45) titleGRAPH <- "Histograma 3D" ################################################################################### library(plot3D) par(mfrow=c(1,1), mar=c(3.5, 3.5, 2, 1), mgp=c(2.25, 0.8, 0)) # # Duplicate gtd -> gtdnew # gtdnew <- data.frame(gtd[,varlistindex]) names(gtdnew)=names(gtd)[varlistindex] # # Calculations and graphics # minX <- min(gtdnew[,1]) maxX <- max(gtdnew[,1]) minX=minX-0.001 maxX=maxX+0.001 minY <- min(gtdnew[,2]) maxY <- max(gtdnew[,2]) minY=minY-0.001 maxY=maxY+0.001 breaksX = c(seq(minX,maxX,by=((maxX-minX)/nclasses[1]))) breaksY = c(seq(minY,maxY,by=((maxY-minY)/nclasses[2]))) breaksX1 = c(seq(minX,maxX,by=((maxX-minX)/(nclasses[1]-1)))) breaksY1 = c(seq(minY,maxY,by=((maxY-minY)/(nclasses[2]-1)))) x_c <- cut(gtdnew[,1], breaks=breaksX) y_c <- cut(gtdnew[,2], breaks=breaksY) # # Calculate joint counts at cut levels # z <- table(x_c, y_c) # # Plot as a 3D histogram # hist3D(x=breaksX1, y=breaksY1, z=z, border="black", theta=angles[1], phi=angles[2], axes=TRUE, facets=TRUE, ticktype="detailed", col = jet.col(100, alpha = 0.75), space=spaceBARS, label=TRUE, nticks=5, scale = T, xlab=paste(names(gtdnew)[1],titleUNITS[1]), ylab=paste(names(gtdnew)[2],titleUNITS[2]), zlab="frequências absolutas", main=titleGRAPH) z