Introduction to Medical Statistics 2026
Exercise class III
Statistical Analysis: Main Concepts and Principles; Binomial Distribution

Author

Ronald Geskus

Published

March 24, 2026

I. Calculation of binomial probabilities

The US CDC estimates that 90% of Americans have had chickenpox by the time they reach adulthood.

Suppose we take a random sample of 100 American adults. Is the use of the binomial distribution appropriate for calculating the probability that exactly 97 out of 100 randomly sampled American adults had chickenpox during childhood? Explain your answer.

Answer: It is. The binomial distribution with “success” probability \(\pi=0.9\) and \(N=100\) describes the variation in number of observed adults that had chickenpox during childhood in a sample of 100.

How many American adults that had chickenpox would you expect to observe among the 100? Wat is the variance of the number of observed cases among 100 sampled individuals? Use the dbinom function to calculate the probability that exactly 97 out of 100 randomly sampled American adults had chickenpox during childhood, i.e. compute P(X=97) if X \(\sim\) B(0.9,100).

Answer: I would expect \(100 \times 0.9=90\) to have had chickenpox; the variance in the observed number is \(N \times p \times (1-p)=100 \times 0.9 \times 0.1=9\).

\(P(X=97)\) is:

dbinom(97, size=100, prob=0.9)

[1] 0.005891602

What is the probability that exactly 3 out of a new sample of 100 American adults have not had chickenpox in their childhood?

Answer: The same probability as in b. Three not having had chickenpox is the same as 97 having had chickenpox.

Plot the probability function of this binomial distribution. See whether the answer from b. corresponds with the value in the plot. Which number has the highest probability and how large is that probability?

x <- 71:100
y <- dbinom(x,100,0.9)
plot(x, y, lwd=2, type="h")
lines(x[90-70], y[90-70], lwd=2, type="h", col="red")

Answer: In the plot, the answer from b. is the height of the bar at x=97. The most likely number is 90, which is number \(90-70=20\) in the generated vector 71:100 that starts at 71. The probability is P(X=90)=0.132

What’s the probability that we observe a number of adults that is at least 5 lower or higher than 90? Hint: use the pbinom function to compute \(P(X \leq 85)\) and \(P(X \geq 95)\). How do you read this number from the plot that you made in b?

pbinom(85, size=100, prob=0.9) + (1-pbinom(94, size=100, prob=0.9) )

[1] 0.1301499

Answer: In the plot, this is the sum of the bar heights \(\leq 85\) and \(\geq 95\).

We actually calculate the p-value of a test of the null hypothesis that p is 0.9 against the two-sided alternative hypothesis that there is a difference if the observed number of infected individuals is 95. Why?

Answer: we calculate the tail probability, in both directions, that the outcome deviates 5 or more from the expected number under the null hypothesis that p=0.9.

II: Inference for a single proportion

A new chemotherapy has entered phase 2 in drug development. A tumour response is defined as a decrease in tumour size by at least 50% within 6 months after therapy initiation. According to experts, the tumour response probability of the drug should be more than 20% in order to proceed to a phase 3 study.

The statistical test of interest is the test for the tumour response probability \(\pi\) with the null hypothesis \(H_0: \pi=0.2\) (or \(H_0: \pi \leq 0.2\)) versus the alternative \(H_A: \pi >0.2\). The phase 2 trial consists of 50 patients. We observe the tumour response for each patient and then perform the statistical test.

Assume that 16/50 patients show a tumour response. Calculate the corresponding one-sided p-value and describe the result and your conclusion in words. Hint: You can either calculate the p-value yourself using the binomial distribution (as in Exercise I.b or I.e) or you can use the prop.test or binom.test function. (Note that these two functions use different methods; only binom.test gives the same answer as the direct calculation.)

# version 1: using the binomial distribution
1-pbinom(15, size=50, prob=0.2)

[1] 0.03080342

# equivalently:
sum(dbinom(16:50, size=50, prob=0.2))

[1] 0.03080342

# version 2: using binom.test
binom.test(16, 50, p = 0.2, alternative = "greater")


    Exact binomial test

data:  16 and 50
number of successes = 16, number of trials = 50, p-value = 0.0308
alternative hypothesis: true probability of success is greater than 0.2
95 percent confidence interval:
 0.2121041 1.0000000
sample estimates:
probability of success 
                  0.32

# version 3: using prop.test
prop.test(16, 50, p = 0.2, alternative = "greater")


    1-sample proportions test with continuity correction

data:  16 out of 50, null probability 0.2
X-squared = 3.7812, df = 1, p-value = 0.02591
alternative hypothesis: true p is greater than 0.2
95 percent confidence interval:
 0.2145141 1.0000000
sample estimates:
   p 
0.32

Answer: We did not specify when there is enough evidence to reject the null hypothesis. Often 0.05 is chosen as cutoff criterion. This cutoff is called significance level. This would suggest to perform a phase 3 trial. Note however that this is a one-sided alternative hypothesis. Then 0.025 is often suggested as significance level. This level would not suggest to perform a phase 3 trial.

Note that the binomial distribution is not symmetric around \(H_0: p=0.2\). Under the null hypothesis the expected number is 10, and \(P(X \geq 16)\) is not equal to \(P(X \leq 4)\), and \(P(X \geq 16)+P(X \leq 4)\) is equal to 0.0492994. This p-value would be just small enough to suggest a phase 3 trial at the \(\alpha=0.05\) significance level. The binom.test function takes this asymmetry into account.

1-pbinom(15, size=50, prob=0.2)+pbinom(4,size=50,prob=0.2)

[1] 0.04929944

## similarly
binom.test(16, 50, p = 0.2, alternative = "two.sided")


    Exact binomial test

data:  16 and 50
number of successes = 16, number of trials = 50, p-value = 0.0493
alternative hypothesis: true probability of success is not equal to 0.2
95 percent confidence interval:
 0.1952042 0.4669938
sample estimates:
probability of success 
                  0.32

The experts expected a tumour response of 40%. Suppose that the true tumour response probability is indeed 0.4. We reject the null hypothesis if we get a one-sided p-value smaller than 0.025. What is the probability that we reject the null hypothesis if we have a sample of 50 patients (not the current sample, but an arbitrary new sample)? This is the basis of the so-called power or sample size calculation. Hint: First show that for an observed number of tumour responses of 16 the test does not reject the null hypothesis, while for 17 it does. Then calculate the probability of \(P(X\geq 17)\) if the alternative is true, i.e. if p=0.4.

 1-pbinom(15, size=50, prob=0.2)

[1] 0.03080342

# p-value > 0.025
 1-pbinom(16, size=50, prob=0.2)

[1] 0.01444166

 # p-value < 0.025
 # hence we reject the null hypothesis if X>16

# calculate probability of X>=17 if the alternative is true (p=0.4)
1-pbinom(16,size=50,prob=0.4)

[1] 0.8439094

Answer: # Bonus: Graph of the distribution assuming that the null hypothesis p=0.2 is true (black bars) and that the alternative p=0.4 is true (red bars).

library(ggplot2)
binom.prob <- data.frame(x=c(0:50-0.1,0:50+0.1),     prob=c(dbinom(0:50,size=50,prob=0.2),dbinom(0:50,size=50,prob=0.4)),
    hypothesis=rep(c("p=0.2 [null]","p=0.4 [alternative]"),c(51,51)))
head(binom.prob)  # plot first 6  rows of data.frame

     x         prob   hypothesis
1 -0.1 1.427248e-05 p=0.2 [null]
2  0.9 1.784060e-04 p=0.2 [null]
3  1.9 1.092737e-03 p=0.2 [null]
4  2.9 4.370946e-03 p=0.2 [null]
5  3.9 1.283965e-02 p=0.2 [null]
6  4.9 2.953120e-02 p=0.2 [null]

tail(binom.prob)  # plot last six rows of data.frame

       x         prob          hypothesis
97  45.1 2.039565e-13 p=0.4 [alternative]
98  46.1 1.477946e-14 p=0.4 [alternative]
99  47.1 8.385509e-16 p=0.4 [alternative]
100 48.1 3.493962e-17 p=0.4 [alternative]
101 49.1 9.507380e-19 p=0.4 [alternative]
102 50.1 1.267651e-20 p=0.4 [alternative]

# we shift the x-values by 0.1 
#in order to prevent the bars in the figure below for p=0.2 and p=0.4 to overlap
ggplot(binom.prob, aes(x,prob,color=hypothesis)) + 
  geom_segment(aes(xend = x, yend = 0)) +
  geom_vline(xintercept=16.5) + theme(legend.position="top")

III: Analysis of a diagnostic test

A diagnostic study evaluated the Platelia NS1 ELISA assay for diagnosis of dengue. The study participants consisted of 853 children admitted to Children’s Hospital #1, Children’s Hospital #2, or HTD from August 2006 to March 2007. Participants were eligible for entry to the study if they had a history of fever of less than seven days and there was a clinical suspicion of dengue. A patient was classified as having dengue by the reference test (gold standard) if there was RT-PCR detection of DENV RNA in plasma, and/or viral culture and/or serological changes in DENV reactive IgM or IgG levels in paired plasma specimens. The data is stored in dengueNS1.csv, the file dengueNS1_description.txt describes the variables in the data set.

Import the dataset and use the summary function for a data summary.

ns <- read.csv("Data/dengueNS1.csv", as.is=FALSE)
summary(ns)

    studyNo       dayOfIllness   dengue            serotype         ns1     
 Min.   :  1.0   Min.   :1.000   no :122   DENV-1      :345   negative:354  
 1st Qu.: 81.0   1st Qu.:3.000   yes:731   DENV-2      :230   positive:499  
 Median :156.0   Median :3.000             DENV-3      : 66                 
 Mean   :157.2   Mean   :3.476             DENV-4      :  7                 
 3rd Qu.:235.0   3rd Qu.:4.000             undetermined:205                 
 Max.   :325.0   Max.   :6.000

Use the table function to create a cross-table of NS1 result and “true” dengue status; you can add column and row totals via the addmargins function. Calculate the prevalence (proportion) of children with dengue.

addmargins(table(ns$ns1,ns$dengue))

          
            no yes Sum
  negative 122 232 354
  positive   0 499 499
  Sum      122 731 853

Answer: Using the prop.test function also provides the 95% confidence interval

# prevalence:
731/853

[1] 0.8569754

prop.test(731,853)


    1-sample proportions test with continuity correction

data:  731 out of 853, null probability 0.5
X-squared = 433.37, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.8312556 0.8794149
sample estimates:
        p 
0.8569754

Calculate the proportion of children with dengue that test positive in the NS1 assay. This quantity is called the sensitivity “SENS” of the NS1 test. Also calculate the proportion of children without dengue that test negative (the specificity “SPEC” of the test). What do you think of the performance of the test?

Answer: Using the prop.test function also provides the 95% confidence interval. SPEC is high, there are no false positive results. However, about 32% of children with dengue are missed with this test.

# sensitivity
499/731

[1] 0.6826265

prop.test(499,731)


    1-sample proportions test with continuity correction

data:  499 out of 731, null probability 0.5
X-squared = 96.793, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.6473053 0.7159998
sample estimates:
        p 
0.6826265

# specificity
122/122

[1] 1

prop.test(122,122)


    1-sample proportions test with continuity correction

data:  122 out of 122, null probability 0.5
X-squared = 120.01, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.9619553 1.0000000
sample estimates:
p 
1

Calculate the proportion of children that have dengue amongst the ones that test positive in the NS1 assay. This is called the positive predictive value (PPV). Also calculate the proportion of children that do not have dengue amongst the ones that test negative in the NS1 assay. Which characteristic of the test is more clinically relevant, SENS and SPEC or PPV and NPV?

Answer: From a clinical perspective PPV and NPV are more relevant. They quantify the probability to have dengue given the test result. SENS and SPEC quantify the test result given the true dengue status. The true status we do not know in practice.

# PPV
499/499

[1] 1

prop.test(499,499)


    1-sample proportions test with continuity correction

data:  499 out of 499, null probability 0.5
X-squared = 497, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.9904764 1.0000000
sample estimates:
p 
1

# NPV
122/354

[1] 0.3446328

prop.test(122,354)


    1-sample proportions test with continuity correction

data:  122 out of 354, null probability 0.5
X-squared = 33.562, df = 1, p-value = 6.902e-09
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.2956821 0.3970138
sample estimates:
        p 
0.3446328

Is there evidence for a difference in sensitivity of the NS1 ELISA assay between DENV serotypes? Use the tbl_summary function from the gtsummary package (No statistical test is required, just an exploratory data analysis.)

library(gtsummary)
ns.truepos <- subset(ns,dengue=="yes") # select only subjects with dengue
tbl_summary(data=ns,by=serotype,include=c(serotype,ns1))

Characteristic	DENV-1 N = 345¹	DENV-2 N = 230¹	DENV-3 N = 66¹	DENV-4 N = 7¹	undetermined N = 205¹
ns1
negative	72 (21%)	96 (42%)	13 (20%)	3 (43%)	170 (83%)
positive	273 (79%)	134 (58%)	53 (80%)	4 (57%)	35 (17%)
¹ n (%)

## a plain table can be obtained via 
table(ns$ns1,ns$serotype)

          
           DENV-1 DENV-2 DENV-3 DENV-4 undetermined
  negative     72     96     13      3          170
  positive    273    134     53      4           35

IV: difference in proportions

The dataset bmData.csv contains selected variables from 300 patients with confirmed bacterial meningitis. They were randomized to either adjunctive dexamethasone therapy or placebo. Import the dataset and use the summary function for a data summary.

bmData <- read.csv("Data/bmData.csv", as.is=FALSE)
summary(bmData)

     patId               group          age        sex          gcs       
 BM 10  :  1   dexamethasone:143   Min.   :15.00   F: 80   Min.   : 3.00  
 BM 100 :  1   placebo      :157   1st Qu.:31.00   M:220   1st Qu.: 9.00  
 BM 101 :  1                       Median :43.00           Median :13.00  
 BM 108 :  1                       Mean   :43.22           Mean   :11.76  
 BM 109 :  1                       3rd Qu.:54.00           3rd Qu.:15.00  
 BM 110 :  1                       Max.   :91.00           Max.   :15.00  
 (Other):294                                                              
    paresis            op.csf         wc.csf         glu.csf         wc.csf.fup     
 Min.   :0.00000   Min.   : 0.0   Min.   :    1   Min.   :  0.10   Min.   :   16.0  
 1st Qu.:0.00000   1st Qu.:15.0   1st Qu.: 1085   1st Qu.: 10.00   1st Qu.:  352.5  
 Median :0.00000   Median :20.0   Median : 3160   Median : 20.00   Median :  902.5  
 Mean   :0.08333   Mean   :22.9   Mean   : 6252   Mean   : 25.03   Mean   : 2883.4  
 3rd Qu.:0.00000   3rd Qu.:30.0   3rd Qu.: 7830   3rd Qu.: 36.00   3rd Qu.: 2375.0  
 Max.   :1.00000   Max.   :55.0   Max.   :64000   Max.   :100.00   Max.   :84000.0  
                   NA's   :63     NA's   :1                        NA's   :22       
   death.6mo     
 Min.   :0.0000  
 1st Qu.:0.0000  
 Median :0.0000  
 Mean   :0.1033  
 3rd Qu.:0.0000  
 Max.   :1.0000

Create a 2x2 table that compares the number of deaths at 6 months between the two treatment arms, using the function xtabs or table. Is there a difference in the survival status at 6 months of follow-up between the two randomized groups?

xtabs(~group+death.6mo, data=bmData)

               death.6mo
group             0   1
  dexamethasone 134   9
  placebo       135  22

Answer: The percentage that died is larger in the placebo group. However, we need to find out whether this could be due to chance. For this, we need to perform a formal test.

Perform a formal test for a difference in survival status. Formulate the null hypothesis and the alternative hypothesis. Use the chisq.test function.

chisq.test(matrix(c(9,22,134,135), nrow=2))


    Pearson's Chi-squared test with Yates' continuity correction

data:  matrix(c(9, 22, 134, 135), nrow = 2)
X-squared = 4.0154, df = 1, p-value = 0.04509

chisq.test(matrix(c(9,22,134,135), nrow=2), correct=FALSE)


    Pearson's Chi-squared test

data:  matrix(c(9, 22, 134, 135), nrow = 2)
X-squared = 4.8125, df = 1, p-value = 0.02825

Answer: The null hypothesis is that dexamethasone does not change the probability of dying within 6 months. One way to formulate this is to say that treatment arm and survival status are independent. If we define \(\pi_D\) as the death risk for those that receive dexamethasone and \(\pi_P\) for those that receive standard of care (placebo arm in this trial), then we can formulate the null hypothesis as \(H_0: \pi_D=\pi_P\).

Use the prop.test function as alternative. Look at the difference in the way the results are shown compared to the chi-squared test.

prop.test(c(22,9), c(157,143))


    2-sample test for equality of proportions with continuity correction

data:  c(22, 9) out of c(157, 143)
X-squared = 4.0154, df = 1, p-value = 0.04509
alternative hypothesis: two.sided
95 percent confidence interval:
 0.003185485 0.151195167
sample estimates:
    prop 1     prop 2 
0.14012739 0.06293706

prop.test(c(22,9), c(157,143), correct=FALSE)


    2-sample test for equality of proportions without continuity correction

data:  c(22, 9) out of c(157, 143)
X-squared = 4.8125, df = 1, p-value = 0.02825
alternative hypothesis: two.sided
95 percent confidence interval:
 0.009866702 0.144513950
sample estimates:
    prop 1     prop 2 
0.14012739 0.06293706

Answer: The prop.test function gives more detailed information. It adds the 95% confidence interval for the difference and is provides the fraction of individuals who died within 6 months in each of the groups. There is moderately strong suggestion that there is a difference in survival status at 6 months. Note that the p-value differs somewhat depending on whether Yates’ continuity correction is used.

We can also use the gtsummary package to obtain all results in tabular format, using the tbl_summary function.

add_p(tbl_summary(data=bmData, by=group, include=c(group,death.6mo), 
      label=list(death.6mo~"died within 6 months")))

Characteristic	dexamethasone N = 143¹	placebo N = 157¹	p-value²
died within 6 months	9 (6.3%)	22 (14%)	0.028
¹ n (%)
² Pearson’s Chi-squared test