modified version of scatter.dapc for plotting points as individual points

rundapc.R

library(adegenet)

data(H3N2)

dapc1 <- dapc(H3N2, pop=H3N2$other$epid, n.pca=30,n.da=6)

custom.cols <- rep(c("red", "blue"), times = c(1000,903))

scatter(dapc1, col = custom.cols)

scatter.dapc.individual.R

scatter.dapc <- function(x, xax = 1, yax = 2, grp = x$grp,
                         col = seasun(length(levels(grp))),
                         pch = 20, bg = "white", solid = .7,
                         scree.da = TRUE, scree.pca = FALSE,
                         posi.da = "bottomright",
                         posi.pca = "bottomleft",
                         bg.inset = "white",
                         ratio.da = .25, ratio.pca = .25,
                         inset.da = 0.02, inset.pca = 0.02, inset.solid = .5,
                         onedim.filled = TRUE, mstree = FALSE,
                         lwd = 1, lty = 1, segcol = "black",
                         legend = FALSE, posi.leg = "topright",
                         cleg = 1, txt.leg = levels(grp),
                         cstar = 1, cellipse = 1.5, axesell = FALSE,
                         label = levels(grp), clabel = 1, xlim = NULL, ylim = NULL,
                         grid = FALSE, addaxes = TRUE, origin = c(0,0),
                         include.origin = TRUE, sub = "", csub = 1, possub = "bottomleft",
                         cgrid = 1, pixmap = NULL, contour = NULL,
                         area = NULL, label.inds = NULL, ...){
    ONEDIM <- xax==yax | ncol(x$ind.coord)==1


    ## recycle color and pch
    pch <- rep(pch, length(levels(grp)))
    col <- transp(col, solid)
    bg.inset <- transp(bg.inset, inset.solid)

    ## handle grp
    if(is.null(grp)){
        grp <- x$grp
    }

    ## handle xax or yax NULL
    if(is.null(xax)||is.null(yax)){
        xax <- 1L
        yax <- ifelse(ncol(x$ind.coord)==1L, 1L, 2L)
        ONEDIM <- TRUE
    }

    ## handle 1 dimensional plot
    if(!ONEDIM){
        ## set par
        opar <- par(mar = par("mar"))
        par(mar = c(0.1, 0.1, 0.1, 0.1), bg=bg)
        on.exit(par(opar))
        axes <- c(xax,yax)

        ## basic empty plot
        s.class(x$ind.coord[,axes], fac = grp, col = col, cpoint = 0,
                cstar = cstar, cellipse = cellipse, axesell = axesell,
                label = label, clabel = clabel, xlim = xlim, ylim = ylim,
                grid = grid, addaxes = addaxes, origin = origin,
                include.origin = include.origin,
                sub = sub, csub = csub, possub = possub,
                cgrid = cgrid, pixmap = pixmap,
                contour = contour, area = area)

        ## add points
        colfac <- pchfac <- grp
        levels(colfac) <- col
        levels(pchfac) <- pch
        colfac <- as.character(colfac)
        pchfac <- as.character(pchfac)
        if(is.numeric(col)) colfac <- as.numeric(colfac)
        if(is.numeric(pch)) pchfac <- as.numeric(pchfac)

        points(x$ind.coord[,xax],
               x$ind.coord[,yax],
               col = col, pch = pchfac, ...)
        s.class(x$ind.coord[,axes], fac = grp, col = col, cpoint = 0,
                add.plot=TRUE, cstar = cstar, cellipse = cellipse,
                axesell = axesell, label = label,clabel = clabel,
                xlim = xlim, ylim = ylim, grid = grid,
                addaxes = addaxes, origin = origin,
                include.origin = include.origin,
                sub = sub, csub = csub, possub = possub,
                cgrid = cgrid, pixmap = pixmap,
                contour = contour, area = area)

        ## Add labels of individuals if specified. Play around with "air" to get
        ## a satisfactory result.

        if (!is.null(label.inds) & is.list(label.inds)) {
            appendList <- function (x, val) {
                                        # recursevly "bind" a list into a longer list,
                                        # from http://stackoverflow.com/a/9519964/322912
                stopifnot(is.list(x), is.list(val))
                xnames <- names(x)
                for (v in names(val)) {
                    x[[v]] <- if (v %in% xnames && is.list(x[[v]]) && is.list(val[[v]]))
                        appendList(x[[v]], val[[v]])
                    else c(x[[v]], val[[v]])
                }
                x
            }

            do.call("orditorp",
                    c(appendList(list(x = x$ind.coord[, c(xax, yax)],
                                      display = "species"),
                                 label.inds)))
        }

        ## add minimum spanning tree if needed
        if(mstree){
            meanposi <- apply(x$tab,2, tapply, grp, mean)
            D <- dist(meanposi)^2
            tre <- ade4::mstree(D)
            x0 <- x$grp.coord[tre[,1], axes[1]]
            y0 <- x$grp.coord[tre[,1], axes[2]]
            x1 <- x$grp.coord[tre[,2], axes[1]]
            y1 <- x$grp.coord[tre[,2], axes[2]]
            segments(x0, y0, x1, y1, lwd = lwd, lty = lty, col = segcol)
        }

    } else {
        ## set screeplot of DA to FALSE (just 1 bar)
        scree.da <- FALSE

        ## get plotted axis
        if(ncol(x$ind.coord)==1) {
            pcLab <- 1
        } else{
            pcLab <- xax
        }
        ## get densities
        ldens <- tapply(x$ind.coord[,pcLab], grp, density)
        allx <- unlist(lapply(ldens, function(e) e$x))
        ally <- unlist(lapply(ldens, function(e) e$y))
        par(bg=bg)
        plot(allx, ally, type = "n",
             xlab = paste("Discriminant function", pcLab),
             ylab = "Density")
        for(i in 1:length(ldens)){
            if(!onedim.filled) {
                lines(ldens[[i]]$x, ldens[[i]]$y, col = col[i], lwd = 2) # add lines
            } else {
                polygon(c(ldens[[i]]$x, rev(ldens[[i]]$x)),
                        c(ldens[[i]]$y, rep(0,length(ldens[[i]]$x))),
                        col = col[i], lwd = 2, border = col[i]) # add lines
            }
            points(x = x$ind.coord[grp==levels(grp)[i], pcLab],
                   y = rep(0, sum(grp==levels(grp)[i])),
                   pch = "|", col = col[i]) # add points for indiv
        }
    }

    ## ADD INSETS ##
    ## group legend
    if(legend){
        ## add a legend
        temp <- list(...)$cex
        if(is.null(temp)) temp <- 1
        if(ONEDIM | temp<0.5 | all(pch == "")) {
            legend(posi.leg, fill = col, legend = txt.leg,
                   cex = cleg, bg = bg.inset)
        } else {
            legend(posi.leg, col = col, legend = txt.leg, cex = cleg,
                   bg = bg.inset, pch = pch, pt.cex = temp)
        }
    }

    ## eigenvalues discriminant analysis
    if(scree.da && ratio.da>.01) {
        inset <- function(){
            myCol <- rep("white", length(x$eig))
            myCol[1:x$n.da] <- "grey"
            myCol[c(xax, yax)] <- "black"
            myCol <- transp(myCol, inset.solid)
            barplot(x$eig, col=myCol, xaxt="n", yaxt="n", ylim=c(0, x$eig[1]*1.1))
            mtext(side=3, "DA eigenvalues", line=-1.2, adj=.8)
            box()
        }

        add.scatter(inset(), posi = posi.da, ratio = ratio.da,
                    bg.col = bg.inset, inset = inset.da)
    }

    ## eigenvalues PCA
    if(scree.pca && !is.null(x$pca.eig) && ratio.pca>.01) {
        inset <- function(){
            temp <- 100* cumsum(x$pca.eig) / sum(x$pca.eig)
            myCol <- rep(c("black","grey"), c(x$n.pca, length(x$pca.eig)))
            myCol <- transp(myCol, inset.solid)
            plot(temp, col=myCol, ylim=c(0,115),
                 type="h", xaxt="n", yaxt="n", xlab="", ylab="", lwd=2)
            mtext(side=3, "PCA eigenvalues", line=-1.2, adj=.1)
        }
        add.scatter(inset(), posi = posi.pca, ratio = ratio.pca,
                    bg.col = bg.inset, inset = inset.pca)
    }


    return(invisible(match.call()))
}