fdiscgammagpd: Fit the discrete gamma-GPD spliced threshold model

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

Description

This function takes count data and fits the gamma-GPD spliced threshold model to it. The model consists of a discrete truncated gamma as the bulk distribution, up to the threshold, and a discrete GPD at and above the threshold. The 'shift' is ideally the minimum count in the sample.

Usage

fdiscgammagpd(x, useq, shift = NULL, pvector=NULL,
std.err = TRUE, method = "Nelder-Mead", ...)

Arguments

x	A vector of count data.
useq	A vector of possible thresholds to search over. These should be discrete numbers.
shift	The amount the distribution is shifted. It is recommended to use the minimum number in the count data when modeling the clone size distribution of the TCR repertoire.
pvector	A vector of 5 elements corresponding to the initial parameter estimates. These 5 initial values are for the gamma shape and rate, the threshold, and the GPD sigma and xi. If they are not prespecified, the function computes pvector automatically.
std.err	Logical. Should the standard errors on the estimates be computed from the Hessian matrix?
method	Character string listing optimization method fed to optim. Defaults to Nelder-Mead.
...	Other arguments passed to the function.

Value

x	Numerical vector of the original data input
shift	Numeric specifying the original shift input.
init	Numerical vector of the initial values of the parameter estimates. This is the same as pvector.
useq	Numerical vector containing the thresholds the grid search was performed over.
nllhuseq	Numerical vector of negative log likelihoods computed at each threshold in useq.
optim	Output from optim for the bulk and tail distributions.
nllh	The negative log likelihood corresponding to the maximum likelihood fitted distribution.
mle	A numerical vector containing the estimates for phi, shape, rate, threshold, sigma, and xi.
fisherInformation	The Fisher information matrix computed from the Hessian output from optim.

Author(s)

[e-mail]

Examples

data("repertoires")
thresholds1 <- unique(round(quantile(repertoires[[1]], c(.75,.8,.85,.9,.95))))
thresholds2 <- unique(round(quantile(repertoires[[2]], c(.75,.8,.85,.9,.95))))

fit1 <- fdiscgammagpd(repertoires[[1]], useq = thresholds1,
                        shift = min(repertoires[[1]]))
fit2 <- fdiscgammagpd(repertoires[[2]], useq = thresholds1,
                        shift = min(repertoires[[2]]))

fit1
fit2

hillarykoch/powerTCR documentation built on March 17, 2021, 8:05 p.m.