R/desponds_functions.R

In powerTCR: Model-Based Comparative Analysis of the TCR Repertoire

# This function does model fitting and threshold selection as described in
# Desponds et al (2016). Details on page 5 of supplement (right column)
# x is a vector of clone sizes

fdesponds <- function(x){
    if(!is(x, "numeric")){
        stop("x must be numeric.")
    }
    
    if(any(x != round(x))){
        stop("all elements in x must be integers.")
    }
    
    x <- sort(x)
    Cmins <- unique(x)
    alpha <- rep(NA, length(Cmins)-1)
    KS <- rep(NA, length(Cmins)-1)

    n <- rep(NA, length(Cmins)-1)
    for(i in seq_along(n)){
        d <- x[x > Cmins[i]]

        n[i] <- length(d)
        alpha[i] <- n[i]*(sum(log(d/Cmins[i])))^(-1) + 1

        empir.cdf <- rep(NA, length(d))
        for(j in seq_along(d)){
            empir.cdf[j] <- mean(d <= d[j])
        }

        est.cdf <- VGAM::ppareto(d, scale=Cmins[i], shape = alpha[i]-1)

        KS[i] <- max(abs(est.cdf - empir.cdf))
    }

    min.KS <- min(KS, na.rm = TRUE)
    Cmin <- Cmins[KS == min(KS, na.rm = TRUE)][1]
    powerlaw.exponent <- alpha[KS == min(KS, na.rm = TRUE)]
    powerlaw.exponent <- powerlaw.exponent[!is.na(powerlaw.exponent)]
    out <- c(min.KS, Cmin, powerlaw.exponent, powerlaw.exponent-1)
    names(out) <- c("min.KS", "Cmin", "powerlaw.exponent", "pareto.alpha")
    out
}