My primary research focuses on developing methodology for statistical inference with complex longitudinal data in comparative effectiveness research. My areas of methodological interest include causal inference, Bayesian statistics, longitudinal data analysis, and mixture and joint modelling.
You can share your program, learning goals, and hobby.
Most of us were trained in the frequentist tradition. Before we introduce a new way of thinking, let us revisit what we already know.
Given the p-value = 0.012 from the binomial exact test, which of the following are true?
Given the 95% CI of \(\theta\) is (0, 0.14), which of the following are true?
The frequentist framework is powerful and well-established. But after seeing that output, does it answer the questions a clinician or clinical researcher actually has?
Consider what we really want to know in the Phase I example:
- What is the probability that the true AE rate is below 20%, given what we just observed?
- What is the probability this drug is safe enough to escalate?
- If I had some prior knowledge about this drug class, how should that change my conclusion?
- A patient asks: “What are the chances this treatment is safe for me?” and what does your analysis let you say?
The honest frequentist answer to all four questions is: I cannot tell you. Not because the methods are wrong, but because they were never designed to answer these questions.
| Frequentist | What we actually want | |
|---|---|---|
| Parameter \(\theta\) | Fixed unknown constant | A quantity we can assign probability to |
| Prior knowledge | Informal; Hypothesis driven | Formally incorporated |
| Direct probability statement about \(\theta\)? | ❌ Not possible, what we get: \(P(\text{data this extreme} \mid H_0 \text{ true})\) | ✅ This is what we want, e.g., \(P( parameter > MCID \mid \text{data})\) |
| 95% CI | 95% of future intervals contain \(\theta\) | 95% probability \(\theta\) is in this interval |
The p-value tells us how surprised we should be by the data if the null were true. It does not tell us how likely the null is to be true.
As Wasserstein and Lazar (2016) note, “p-values do not measure the probability that the studied hypothesis is true.”
The Bayesian framework is built around a simple and intuitive idea:
We start with a belief. We observe data. We update our belief.
\[\text{updated belief} = \text{current evidence} \times \text{prior belief or evidence}\]
This is not a new idea, it is how humans naturally reason under uncertainty.
Bayesian statistics formalizes this process using probability.
Returning to our Phase I example:
Before seeing data: we have some belief about \(\theta\), the true AE rate. Perhaps we think all values between 0 and 1 are equally plausible (a flat, uninformative prior). Or perhaps we have some prior clinical experience suggesting the rate is likely low.
We observe data: 0 out of 20 patients experienced an AE.
We update our belief: combining what we knew before with what we just observed, we obtain a posterior distribution for \(\theta\). This is a full probability distribution over all possible values of \(\theta\), not just a point estimate.
We answer the clinical question directly: from the posterior, we can compute \(P(\theta < 0.20 \mid \text{data})\) directly. This is the probability that the true AE rate is below 20%, given everything we know.
More about Beta-binomial model next week. For illustration, from the posterior distribution \(Beta(1, 21)\), we can directly compute:
# Posterior: Beta(1, 21) given prior Beta(1,1) and x=0, n=20
a_post <- 1 + 0
b_post <- 1 + 20
# Posterior mean
cat("Posterior mean:", round(a_post / (a_post + b_post), 3), "\n")Posterior mean: 0.045 # 95% Credible Interval
cat("95% Credible Interval:",
round(qbeta(c(0.025, 0.975), a_post, b_post), 3), "\n")95% Credible Interval: 0.001 0.161 # Direct answer to the clinical question
cat("P(theta < 0.20 | data):",
round(pbeta(0.20, a_post, b_post), 3), "\n")P(theta < 0.20 | data): 0.991 \(P(\theta < 0.20 \mid \text{data}) = 0.99\) means there is a 99% posterior probability that the true AE rate is below 20%. We can now directly and confidently recommend dose escalation.
This is the kind of direct probability statement that clinical researchers and regulators find intuitive and actionable.
Bayesian methods are not new, they date back to the Reverend Thomas Bayes in the 18th century. But their adoption in clinical research has accelerated dramatically in recent decades, driven by three forces:
1. Computational advances. Modern Markov Chain Monte Carlo (MCMC) methods and software (Stan, JAGS, brms in R) have made Bayesian analysis of complex models practical for applied researchers.
2. Regulatory acceptance. The FDA Center for Devices and Radiological Health finalized Bayesian guidance for medical device trials in 2010. Bayesian adaptive designs have been increasingly used in oncology, rare disease, and pediatric trials where small sample sizes make frequentist inference particularly limited (Lee et al. 2026).
3. Recognition of frequentist limitations. The American Statistical Association’s 2016 statement on p-values (Wasserstein and Lazar 2016) and its 2019 follow-up (Wasserstein, Schirm, and Lazar 2019) sparked widespread reflection on over-reliance on p-values and arbitrary significance thresholds, opening the door for Bayesian alternatives.
By the end of this course, you will have the tools to read, critically evaluate, and conduct Bayesian analyses in your own clinical research.
1. Long run relative frequency
The classical or frequentist way of interpreting probability is based on the proportion of times an outcome occurs in a large (i.e. approaching infinity) number of repetitions of an experiment.
2. Subjective probability
The alternative way (essentially proposed by Bayes, LaPlace, other later prominent Bayesians) is the subjective probability.
Probability has a central place in Bayesian analysis
Sample space This set of all possible outcomes of an experiment/trial is known as the sample space of the experiment/trial and is denoted by \(S\) or \(\Omega\).
Experiment/Trial: each occasion we observe a random phenomenon that we know what outcomes can occur, but we do not know which outcome will occur
Event: Any subset of the sample space \(S\) is known as an event, denoted by \(A\). Note that, \(A\) is also a collection of one or more outcomes.
Probability defined on events: For each event \(A\) of the sample space \(S\), \(P(A)\), the probability of the event \(A\), satisfies the following three conditions:
Union: Denote the event that either \(A\) or \(B\) occurs as \(A\cup B\).
Intersection: Denote the event that both \(A\) and \(B\) occur as \(A\cap B\)
Complement: Denote the event that \(A\) does not occur as \(\bar{A}\) or \(A^{C}\) or \(A^\prime\) (different people use different notations)
Disjoint (or mutually exclusive): Two events \(A\) and \(B\) are said to be disjoint if the occurrence of one event precludes the occurrence of the other. If two events are mutually exclusive, then \(P(A\cup B)=P(A)+P(B)\)
Sample Space
Event
Probability of an event. Let \(R\) be the sum of two standard dice. Suppose we are interested in \(P(R \le 4)\). Notice that the pair of dice can fall 36 different ways (6 ways for the first die and six for the second results in 6x6 possible outcomes, and each way has equal probability \(1/36\). \[\begin{aligned} P(R \le 4 ) &= P(R=2)+P(R=3)+P(R=4) \\ &= P(\left\{ 1,1\right\} )+P(\left\{ 1,2\right\} \mathrm{\textrm{ or }}\left\{ 2,1\right\} )+P(\{1,3\}\textrm{ or }\{2,2\}\textrm{ or }\{3,1\}) \\ &= \frac{1}{36}+\frac{2}{36}+\frac{3}{36} \\ &= \frac{6}{36} \\ &= \frac{1}{6} \end{aligned}\]
General Addition Rule: \(P(A\cup B)=P(A)+P(B)-P(A\cap B)\)
The reason behind this fact is that if there is if \(A\) and \(B\) are not disjoint, then some area is added twice when we calculate \(P\left(A\right)+P\left(B\right)\). To account for this, we simply subtract off the area that was double counted. If \(A\) and \(B\) are disjoint, \(P(A\cup B)=P(A)+P(B)\).
Complement Rule: \(P(A)+P(A^c)=1\)
This rule follows from the partitioning of the set of all events (\(S\)) into two disjoint sets, \(A\) and \(A^c\). We learned above that \(A \cup A^c = S\) and that \(P(S) = 1\). Combining those statements, we obtain the complement rule.
Completeness Rule: \(P(A)=P(A\cap B)+P(A\cap B^c)\)
This identity is just breaking the event \(A\) into two disjoint pieces.
Law of total probability (unconditioned version): Let \(A_1, A_2, \ldots\) be events that form a partition of the sample space \(S\), Let \(B\) be any event. Then, \[P(B) = P(A_1 \cap B) + P(A_2 \cap B) + \ldots. \]
This law is key in deriving the marginal event probability in Bayes’ rule. Recall the HIV example from session 1, we have \(P(T^+) = P(T^+ \cap D^+)+P(T^+ \cap D^-)\).
Conditional Probability: The probability of even \(A\) occurring under the restriction that \(B\) is true is called the conditional probability of \(A\) given \(B\). Denoted as \(P(A|B)\).
In general we define conditional probability (assuming \(P(B) \ne 0\)) as \[P(A|B)=\frac{P(A\cap B)}{P(B)}\] which can also be rearranged to show \[\begin{aligned} P(A\cap B) &= P(A\,|\,B)\,P(B) \\ &= P(B\,|\,A)\,P(A) \end{aligned}\]
Suppose I know that whenever it rains there is 10% chance of being late at work, while the chance is only 1% if it does not rain. Suppose that there is 30% chance of raining tomorrow; what is the chance of being late at work?
Denote with event \(A\) as “I will be late to work tomorrow” and event \(B\) as “It is going to rain tomorrow”
\[\begin{aligned} P(A) &= P(A\mid B)P(B) + P(A \mid B^c) P(B^c) \\ &= 0.1\times 0.3+0.01\times 0.7=0.037 \end{aligned}\]
Independent: Two events \(A\) and \(B\) are said to be independent if \(P(A\cap B)=P(A)P(B)\).
What independence is saying that knowing the outcome of event \(A\) doesn’t give you any information about the outcome of event \(B\). If \(A\) and \(B\) are independent events, then \(P(A|B) = P(A)\) and \(P(B|A) = P(B)\).
In simple random sampling, we assume that any two samples are independent. In cluster sampling, we assume that samples within a cluster are not independent, but clusters are independent of each other.
Assumptions of independence and non-independence in statistical modelling are important.
In linear regression, for example, correct standard error estimates rely on independence amongst observations.
In analysis of clustered data, non-independence means that standard regression techniques are problematic.
Bayes’ Rule: This arises naturally from the rule on conditional probability. Since the order does not matter in \(A \cap B\), we can rewrite the equation:
\[\begin{aligned} P(A\mid B) &= \frac{P(A \cap B) }{P(B)}\\ &= \frac{P(B \cap A)}{P(B)}\\ &= \frac{P(B\mid A)P(A) }{P(B)}\\ &= \frac{P(B\mid A)P(A) }{P(B \cap A) + P(B \cap A^C)}\\ &= \frac{P(B\mid A)P(A) }{P(B\mid A)P(A) + P(B\mid A^C)P(A^C)} \end{aligned}\]
Frequency-type (Empirical Probability): based on the idea of frequency or long-run frequency of an event.
\[ P(A) = \frac{|A|}{|S|} = \frac{\text{Number of outcomes in S that satisfy A}}{\text{Total number of outcomes in S} } \]
n = 1000
pHeads =0.5
set.seed(123)
flipSequence = sample( x=c(0,1),prob = c(1-pHeads,pHeads),size=n,replace=T)
r = cumsum(flipSequence)
n= 1:n
runProp = r/n
flip_data <- data.frame(run=1:1000,prop=runProp)
ggplot(flip_data,aes(x=run,y=prop,frame=run)) +
geom_path()+
xlim(1,1000)+ylim(0.0,1.0)+
geom_hline(yintercept = 0.5)+
ylab("Proportion of Heads")+
xlab("Number of Flips")+
theme_bw()
Belief-type (Subjective Probability): based on the idea of degree of belief (weight of evidence) about an event where the scale is anchored at certain (=1) and impossible (=0).
Which probability to use?
For inference, the Bayesian approach relies mainly on the belief interpretation of probability.
The laws of probability can be derived from some basic axioms that do not rely on long-run frequencies.
Random variables are used to describe events and their probability. The formal definition just means that a random variable describes how numbers are attached to each possible elementary outcome. Random variables can be divided into two types:
Continuous Random Variables: the variable takes on numerical values and could, in principle, take any of an uncountable number of values. These typically have values that cannot be counted, e.g. age, height, blood pressure, etc.
Discrete Random Variables: the variable takes on one of small set of values (or only a countable number of outcomes).
The value of a particular measurement can be thought of as being determined by chance, even though there may be no true randomness
Probability distributions
“All models are wrong but some are useful.”
This is a function with two important properties
The heights of the function at two different x-values indicate relatively how likely those x-values are.
The area under the curve between any two x-values is the probability that a randomly sampled value will fall in that range
The area under the whole curve is equal to the total probability, i.e. equal to 1
From histograms to PDFs
With a large enough sample, and narrow enough intervals, histogram starts to look smooth.
The relative heights of the histogram at each blood pressure value tell us how likely these values are.
For example, near 120 mm Hg there are about 800 people in each 5 mm interval.
As the width of the interval approaches zero, the probability of being in that interval approaches zero. e.g. The event that “BP is exactly 129” has zero probability, with \[P(BP = 129) \approx P(128.99999 < BP < 129.00001) \approx 0.\]
For a continuous random variable, we describe the shape of the distribution with a smooth curve (pdf) instead of a histogram.
The relative height of the curve at any two points indicates the relative likelihood of the two values. The height is not a probability.
The area under the curve between any two points is the probability that a randomly sampled value will fall in that range. e.g., the shaded area represent = 0.5. In fact, the probability of a BP from this population being between 112 and 145 is 0.5.
Discrete random variables and distributions means that the random variable can take on only a countable number of values. We will describe these distributions using probability mass function.
Continuous random variables and distributions: We will describe these distributions using probability density function.
The expectation for discrete distributions \[E[X]= \sum_{x=0}^{n}x\,P(X=x) \]
and the variance for discrete distributions \[V[X]=\sum_{x=0}^{n}\left(x-E\left[X\right]\right)^{2}\,P(X=x) \]
The expectation for continuous distributions \[E[X]= \int_{a}^{b}x\,f(x) dx\]
and the variance for continuous distributions \[V[X]= E[(X-E(X))^2] = E(X^2) - E(X)^2\] and \[E[X^2]= \int_{a}^{b}x^2\,f(x) dx\]
Daniel Bernoulli FRS (1700-1782) was a Swiss mathematician and physicist and was one of the many prominent mathematicians in the Bernoulli family from Basel.
A binomial experiment is one that:
The Binomial Probability Mass Function
The probability mass function for the binomial distribution is \[ P(Y=y) = \binom{n}{y} p^y(1-p)^{n-y}, \ \text{for } y=0,\ldots,n \] where \(n\) is a positive integer and \(p\) is a probability between 0 and 1.
This probability mass function uses the expression \(\binom{n}{y}\), called a binomial coefficient \[ \binom{n}{y} = \frac{n!}{y!(n-y)!} \] which counts the number of ways to choose \(x\) things from \(n\).
The mean of the binomial distribution is \(E(Y) = np\) and the variance is \(Var(Y) = np(1-p)\).
We can use R to calculate binomial probabilities. In general, for any distribution, the “d-function” gives the density function (pmf or pdf). \(P( x = 5 | n=10, p=0.8 ) = 0.02642412\) from dbinom(5, size=10, prob=0.8) and \(P( x = 8 | n=10, p=0.8 )=0.3019899\) from dbinom(8, size=10, prob=0.8).
A commonly used distribution for count data is the Poisson. e.g., number of patients arriving at ER over a 12 hr interval.
The following conditions apply:
Assuming that these conditions hold for some count variable \(Y\), the the probability mass function is given by \[P(Y=y)=\frac{\lambda^{y}e^{-\lambda}}{y!}, y = 0, 1, \ldots\] where \(\lambda\) is the expected number of events over 1 unit of time or space and \(e\) is the euler’s number \(2.718281828\dots\).
\[E[Y] = \lambda\] \[Var[Y] = \lambda\]
Example: Suppose we are interested in the population size of small mammals in a region. Let \(Y\) be the number of small mammals caught in a large trap over a 12 hour period. Finally, suppose that \(Y\sim Poisson(\lambda=2.3)\). What is the probability of finding exactly 4 critters in our trap? \[P(Y=4) = \frac{2.3^{4}\,e^{-2.3}}{4!} = 0.1169\] What about the probability of finding at most 4? \[\begin{aligned} P(Y\le4) &= P(Y=0)+P(Y=1)+P(Y=2)+P(Y=3)+P(Y=4) \\ &= 0.1003+0.2306+0.2652+0.2033+0.1169 \\ &= 0.9163 \end{aligned}\]
What about the probability of finding 5 or more? \[P(Y\ge5) = 1-P(Y\le4) = 1-0.9163 = 0.0837\]
These calculations can be done using the density function (d-function) for the Poisson and the cumulative distribution function (p-function). e.g. \(P(Y=4\mid \lambda = 2.3\) can be obtained using dpois(4, lambda=2.3) and \(P(Y \leq 4\mid \lambda = 2.3\) can be obtained using ppois(4, lambda=2.3).
If \(Y\) is the count of independent Bernoulli trials required to achieve the \(r\)th successful trial (\(k\) failures) when the probability of success is \(p\).
\(Y \sim NB(r, p)\) and the probability mass function is
\[ f(Y = k) = \frac{\Gamma(k+r)}{k! \Gamma(r)} (1-p)^y (p)^{r}, \ \text{for } y=0,1,\ldots \]
\(E(Y) = \frac{pr}{1-p}\) and \(V(Y) = \frac{pr}{(1-p)^2}\)
Connection between Poisson and Negative Binomial.
Let \(\lambda = E(Y) = \frac{pr}{1-p}\), we thus have \(p=\frac{\lambda}{r+\lambda}\) and \(V(Y) = \lambda(1 + \frac{\lambda}{r}) > \lambda\).
The variance of the NB is greater than its mean, where as the variance of Poisson equals to its mean!
When we fit a Poisson regression and encounter over-dispersion - the observed variance of the data is greater than the variance of the Poisson model, we can fit a Negative Binomial regression instead!
Suppose you wish to draw a random number number between 0 and 1 and any two intervals of equal size should have the same probability of the value being in them. This random variable is said to have a Uniform(0,1) distribution.
Undoubtedly the most important distribution in statistics is the normal distribution. Normal Distributions has two parameters: mean \(\mu\) and standard deviation \(\sigma\). Its probability density function is given by \[f(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right]\] where \(\exp[y]\) is the exponential function \(e^{y}\). We could slightly rearrange the function to
\[f(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^{2}\right]\]
and see this distribution is defined by its expectation \(E[X]=\mu\) and its variance \(Var[X]=\sigma^{2}\).
Example: It is known that the heights of adult males in the US is approximately normal with a mean of 5 feet 10 inches (\(\mu=70\) inches) and a standard deviation of \(\sigma=3\) inches. I am 5 feet 4 inches (64 inches). What proportion of the adult male population is shorter than me?
Using R you can easily find this probability is 0.02275013 using pnorm(64, mean=70, sd=3).
The probability density function \[ f(x) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} x^{\alpha - 1} (1-x)^{\beta - 1}, 0 \leq x \leq 1.\] where \[ B(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}\] and \(\Gamma()\) is called the gamma function.
The probability density function \[f(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x},x \geq 0.\]
Role of these distributions
All of them are potentially distributions for observed outcomes. Similar choices as in non-Bayesian modelling
A select group of continuous distributions are frequently used to represent prior distributions for parameters