Van Andel Institute
Suppose I want to model whether something happens (1 if so, 0 otherwise). For example, whether someone is admitted to a UC Berkeley graduate program. We might want to model the importance of factors to the outcome probability \(p\). One way to do this is with a binomial variable \(X \sim Bin(n, p)\), where the number of individuals applying is \(n\), the probability of admission is \(p\), and \(X \sim Bin(n, p)\), where \(p\) might depend on gender, department & skills. Skills for graduate research are notoriously difficult to quantify, so we might instead ask whether gender bias is uncomfortably good at predicting admission.
\[ \begin{align} Pr(X = k | n, p) &= \binom{n}{k} p^{k} \left( 1 - p \right)^{n-k} \\ \binom{n}{k} &= \frac{n!}{k! \left( n - k \right)!} \\ E\left[X\right] &= np \\ Var\left[X\right] &= np(1-p) \end{align} \]
Recall from the previous slide that a series of successes or failures can be modeled as \(X \sim Bin(n, p)\). If we model \(p\) and use the logistic function, \(\text{logit}(p) = log\frac{p}{1-p}\), we can perform logistic regression.
Since \(\text{logit}(p)\) can potentially be infinite, we can’t perform a least- squares fit the way we might for a linear model. But we can finesse this by maximizing the likelihood (joint probability of seeing the observed data). In the next slide we will see an example. For more on maximum likelihood and logistic models, see this notebook. For much more on both, see marginaleffects.com.
UCBAdmissions is included with R, in the datasets library. Data is provided as frequencies, so we use weights in glm().
It turns out that many (perhaps most) categorical outcomes can be modeled with a flexible, linear-like model. This is particularly handy if predictors interact in complicated ways. To implement it, we must predict how many events occur.
A rule of thumb is that if \(X ~ Bin(n, p)\), \(n \geq 20\), and \(p \leq .05\), then \(X \approx Poisson(\lambda = np)\). This is handy for modeling counts of events. Let’s see how moving from logistic to Poisson models could help us do so.
Suppose we have events that occur at varying times but at a steady rate.
A classic example is radioactive decay. At any instant, I can’t tell you if a particle will be detected. But if you tell me the isotopic half-life, I can guess how many events you’ll see in a day, year, or century. If the intensity of decay is \(\lambda\), I can make good guesses about their counts \(X\) over time.
We say \(X \sim Pois(\lambda)\): “\(X\) is Poisson with intensity \(\lambda\)”.
\[ \begin{align} E\left[X\right] &= \lambda \\ Var\left[X\right] &= \lambda \\ Pr(X=k|\lambda) &= \frac{\lambda^{k}}{k!}e^{-\lambda} \end{align} \]
Suppose instead of number of events, \(X\), we track time between events, \(Y\).
Then \(Y \sim Exp(\lambda)\): “\(Y\) is exponential with intensity \(\lambda\)”.
\[ E\left[Y\right] = \frac{1}{\lambda} \]
The probability of \(k\) events per observation \(X\) is \(Pr(\sum^{k}_{i=1}Y_i \leq 1)\):
\[ \begin{align} Pr(0|\lambda) &= e^{-\lambda} \\ Pr(1|\lambda) &= \lambda e^{-\lambda} \\ Pr(2|\lambda) &= \frac{\lambda^2}{2!} e^{-\lambda} \\ ... \\ Pr(k|\lambda) &= \frac{\lambda^k}{k!} e^{-\lambda} \end{align} \]
If the probability of an event per unit observation \(\lambda = \frac{1}{10}\), the mean number of observations between events will be \(\frac{1}{\lambda} = 10\).
Suppose we have \(j\) exponentially distributed processes \(Y_j \sim Exp(\lambda)\).
If \(Z = \sum^{j}_{i=1} Y_j\), we can model the sum as \(Z \sim Gamma(j, \lambda)\):
\[ \begin{align} E\left[Z\right] &= \sum^{j}_{i=1} E\left[Y_j\right] = \frac{j}{\lambda} \\ Var\left[Z\right] &= \sum^{j}_{i=1} Var\left[Y_j\right] = \frac{j}{\lambda^2}\\ Pr(Z=z|j,\lambda) &= \frac{z^{j-1}e^{-z/\lambda}}{\lambda^j\Gamma(j)}, \text{where}\ \Gamma(n) = (n - 1)! \end{align} \]
Recall: if \(X \sim Poisson(\lambda)\), \(t(x_{i+1}) - t(x_i) \sim Exp(\lambda)\).
Now we have the pieces to reconstruct a Poisson process across discrete samples. We can reason about how precise our estimate of \(\lambda\) is, within and across samples, using a model that captures the sampling distribution of \(\lambda\).
The likelihood \(\mathcal{L}(\theta|X)\) of model parameters \(\theta\) is the joint probability of seeing the data \(X\) if \(\theta\) describes the generating process. Note: it’s easier to compute \(\ell(\theta|X) = log(\mathcal{L}(\theta|X))\).
So if \(X \sim Poisson(\lambda)\), the only parameter in \(\theta\) is \(\lambda\). If we label each of \(n\) observations as \(x_i\), then
\[ \begin{align} Pr(x_i|\lambda) &= \frac{\lambda^{x_i}e^{\lambda}}{x_i!}, for i \in 1,\cdots,n\\ \mathcal{L}(\lambda|X) &= Pr(X|\lambda) \\ &= \prod_{i=1}^{n}\frac{\lambda^{x_i}e^{\lambda}}{x_i!}\\ \ell(\lambda|X) &= log(\prod_{i=1}^{n}\frac{\lambda^{x_i}e^{\lambda}}{x_i!}) \\ &= \sum_{i=1}^{n} log( \end{align} \]
It would be nice to have some way to relate this to our sample size \(n\). Luckily for us, we will accomplish that using the Gamma distribution in the next slide.
Suppose that the sole parameter in our Poisson model, \(\lambda\), itself has some randomness to it. A natural model for this is the Gamma distribution:
\[ \begin{align} \end{align} \]
It would be nice to have some way to relate this to our sample size \(n\). Luckily for us, we will accomplish that using the Gamma distribution in the next slide.
Suppose we look at inherited genetic variants in the general population; in a disease; and in a specific type or age group of disease. We might wonder whether the rates of inherited genetic variants (of any sort!) are in each group, and we might wonder how strong the evidence is for differences between groups. But if the disease is rare, we know that our estimate of \(\lambda\) cannot be as precise as it would be if the disease was common, because we simply don’t have as many observations. The Gamma-Poisson distribution helps us represent this inherent imprecision in our simulations and reason about how that affects our comparisons.
Let’s make our examples concrete by looking at germline genetic variants in an unselected control population, versus genetic variants in leukemia patients. We will assume that clinically relevant genetic variants are somewhat rare, but previous studies have reported that at least 5% of leukemia patients harbor inherited genetic variants which are rare in the general population.
Recall that one type of Gamma distribution is the sum of exponential processes.
A statistical model for graph traversal is the exponential random graph model, or Erdos-Renyi model, where the number of edges (connections) from any given node (such as a gene, or a protein-protein interaction) is exponential.
Since the sum of traversal times for disconnected graph components is a sum of exponential processes under this random graph model, we can test the hypothesis that the intensity of emissions (“hops” before hitting a gene variant, whether inherited or acquired) is the same across disease subtypes, and we can do this with an amended graph that reflects any protein-protein fusions in a patient:
(overdispersed loglinear models)