The best way to plot high amount of discrete data with 2 variables in R I am trying to see the relationship between two variables (say A and B) in a plot in RStudio. Both are discrete and range from 1 to 10. However, I have 1000s of data points, so given that there are only 100 possible spaces in which there can be a point, almost every possible place on a graph has a point.
How can I represent 1000s of points on a plot like this, whilst being able to see how many are at each point too?
 A: Mosaic plots are a good way of doing this.
https://cran.r-project.org/web/packages/ggmosaic/vignettes/ggmosaic.html
A: The ggplot2 library should handle something like this. There are example of the specific code out on the internet. I’ll just address the idea, since this is CV.SE, not SO.
I would represent the points in a data frame with three columns. One column would have the x-coordinate, one column would have the y-coordinate, and one column would have the count of how many instances of that x-y pair there are. Then you can use a color to denote the prevalence of a point, which ggplot2 can do.
A: One potential option is to add a tiny bit of random noise to each observation. In that way fewer points will overlap.
You can either add it directly and use R's basic plotting capabilities or look into the jitter type layer that comes with the GGplot package that adds the noise automatically.
A: Similar to what Dave proposes, but in base R: visualize table counts using grayscale, with darker grays for cells with higher counts.

set.seed(1)
nn <- 1e6
aa <- sample(1:10,nn,prob=(1:10)^2-5*(1:10)+20,replace=TRUE)
bb <- sample(1:10,nn,prob=20-(1:10),replace=TRUE)

data_table <- table(aa,bb)

grayscale <- function ( cnt ) paste0("grey",100-3*round(cnt/1000,0))
# this relies on the fact that counts are between 3000 and 30000
# adapt as needed

plot(c(0,12),c(0,11),type="n",las=1,xlab="A",ylab="B")
for ( ii in rownames(data_table) ) {
    for ( jj in colnames(data_table) ) {
        rect(as.numeric(ii)-.5,as.numeric(jj)-.5,as.numeric(ii)+.5,as.numeric(jj)+.5,
            border=NA,col=grayscale(data_table[ii,jj]))
        # optionally, add counts
        # text(as.numeric(ii),as.numeric(jj),data_table[ii,jj],
        #   col=if(data_table[ii,jj]>quantile(data_table,0.7)) "white" else "black")
    }
}
counts_for_legend <- round(seq(min(data_table),max(data_table),length.out=5),0)
legend("right",pch=22,pt.bg=grayscale(counts_for_legend),legend=counts_for_legend,pt.cex=1.5)

Of course, this could be prettified a lot, especially the legend - the question is whether you want to do this by hand (if you want to create this plot only a single time), or programmatically (if this needs to be created often, with different datasets).
Alternatively, if you want a little more color in your life, you could change the grayscale() function above to one that outputs a black body radiation color:

lackBodyRadiationColors <- function(x, max_value=1) {
    # x should be between 0 (black) and 1 (white)
    # if large x come out too bright, constrain the bright end of the palette
    #     by setting max_value lower than 1
    foo <- colorRamp(c(rgb(0,0,0),rgb(1,0,0),rgb(1,1,0),rgb(1,1,1)))(x*max_value)/255
    apply(foo,1,function(bar)rgb(bar[1],bar[2],bar[3]))
}

plot(c(0,12),c(0,11),type="n",las=1,xlab="A",ylab="B")
for ( ii in rownames(data_table) ) {
    for ( jj in colnames(data_table) ) {
        rect(as.numeric(ii)-.5,as.numeric(jj)-.5,as.numeric(ii)+.5,as.numeric(jj)+.5,
            border=NA,col=blackBodyRadiationColors(1-data_table[ii,jj]/max(data_table)))
        # optionally, add counts
        # text(as.numeric(ii),as.numeric(jj),data_table[ii,jj],
        #   col=if(data_table[ii,jj]>quantile(data_table,0.7)) "white" else "black")
    }
}
counts_for_legend <- round(seq(min(data_table),max(data_table),length.out=5),0)
legend("right",pch=22,pt.bg=blackBodyRadiationColors(1-counts_for_legend/max(data_table)),
  legend=counts_for_legend,pt.cex=1.5)

