Feb. 26, 2026
2004: BS Math, Kansas State
2006: MS Statistics, Kansas State
2008: Married to Kate Allred
2010: PhD Statistics, University of South Carolina
2010: Assist. Professor, OSU Department of Statistics
2015: Max Habiger born
2021: Professor, Department of Statistics
Other things
3 cats (Marvin is my favorite) and a dog
Basketball, Football, Tennis, Discgolf, Guitar/Piano
From the Abstract: “The rhizosphere is populated by a numerous and diverse array of rhizobacteria, and many impact productivity in largely unknown ways… Operational taxonomic units (OTUs) that were significantly associated with biomass productivity were identified using an exact test adjusted for the false discovery rate…”
url = "https://raw.githubusercontent.com/Jhabige/Jhabige.github.io/refs/heads/master/assets/bacteria.txt" bacteria=read.table(url, row.names=1) biomass=c(0.86, 1.34,1.81,2.37,3.0) colnames(bacteria)= biomass bacteria[1:5,]
## 0.86 1.34 1.81 2.37 3 ## 1847 0 0 4 11 12 ## 588 0 0 0 0 14 ## 1697 0 0 4 1 13 ## 1573 0 0 2 3 7 ## 385 2 0 2 2 12
dim(bacteria)
## [1] 778 5
Shoot biomass \(\approx\) productivity measured for each sample
p = 778 bacteria abundances measured per sample
n = 5 samples
p>n means HIGH DIMENSIONAL data
boxplot(bacteria)
l.bacteria=log(bacteria+1) boxplot(l.bacteria)
biplot(princomp(l.bacteria), scale=0, ylim = c(-3,4), xlim = c(-3,10)) abline(h=0,v=0)
Observe bacteria 508 (top) vs 42, 43 (right). Can we interpret?
What are the arrows?
round(princomp(l.bacteria)$loadings[,1:2],2)
## Comp.1 Comp.2 ## 0.86 0.50 0.85 ## 1.34 0.43 -0.09 ## 1.81 0.47 -0.28 ## 2.37 0.46 -0.39 ## 3 0.37 -0.19
Points: For each bacteria plug in its z score for X1, z score for X2, …
Z1 = .50 X1 + .43 X2 + .47 X3 + .46 X4 + .37 X5 (abundance component)
Z2 = .85 X1 - .09 X2 -.28 X3 -.39 X4 - .19 X5 (association component)
Can you interpret bacteria 508, 42, 43 on the biplot now?
Arrows: Plots of loadings above
interesting=c("508","42","43")
matplot(biomass,t(l.bacteria[interesting,]), pch = 2:4,col=2:4,type="b",ylab="log abundance")
legend("topright",pch=2:4, col=2:4,legend=interesting)
lines(biomass, colMeans(l.bacteria))
Indeed 42, 43 are abundant
Indeed 508 is negatively associated with biomass
plot(princomp(l.bacteria), main = "Variance of PC's")
Abundance variation (Z1) and Productivity associated variation (Z2)
Comp.3 - Comp.5 variation not on biplot!
In Principle Component Analysis:
Z1 and Z2 are linear combinations of X1,..,X5 that maximize variance
Typically PCA is ran after centering and scaling variables to have mean 0 and standard deviation 1.
This way unit of measure doesn’t matter
Not necessary if variables on same scale
Important: Here we ran PCA on samples
Are samples already on “the same scale”, i.e. “technical variation”?
If yes, PCA results won’t change after scaling
boxplot(l.bacteria)
sl.bacteria = scale(l.bacteria) boxplot(sl.bacteria)
round(princomp(l.bacteria)$loadings[,1:2],2)
## Comp.1 Comp.2 ## 0.86 0.50 0.85 ## 1.34 0.43 -0.09 ## 1.81 0.47 -0.28 ## 2.37 0.46 -0.39 ## 3 0.37 -0.19
round(princomp(sl.bacteria)$loadings[,1:2],2)
## Comp.1 Comp.2 ## 0.86 0.37 0.87 ## 1.34 0.46 0.13 ## 1.81 0.50 -0.22 ## 2.37 0.48 -0.36 ## 3 0.41 -0.23
biplot(princomp(l.bacteria), scale=0)
biplot(princomp(sl.bacteria), scale = 0)
matplot(biomass,t(l.bacteria[interesting,]), pch =2:4,col=2:4,type="b",ylab="log abundance")
legend("topright",pch=2:4, col=2:4,legend=interesting)
lines(biomass, colMeans(l.bacteria))
matplot(biomass,t(sl.bacteria[interesting,]), pch = 2:4,col=2:4,type="b",ylab="log abundance")
legend("topright",pch=2:4, col=2:4,legend=interesting)
lines(biomass, colMeans(sl.bacteria))
matplot(biomass,t(sl.bacteria),type="l",ylab="log abundance")
legend("topright",legend="Ummm? Not it!")
The False Discovery Rate (FDR) is \[FDR = E\left[\frac{\# \text{ False Discoveries}}{\#\text{ Discoveries} }\right]\]
Specify statistical model and null hypotheses
Compute p-values (classical approach) \(p_1,p_2,..,p_M\) or posterior null probabilities (empirical Bayes approach) \(q_1,q_2,...,q_M\)
Use statistics above to estimate FDR or “adjust” for FDR
Classical (BH procedure) estimate: \(\widehat{FDR}(p_{(k)}) = \frac{M p_{(k)}}{k}\)
Empirical Bayes estimate: \(\widehat{FDR}(q_{(k)}) = \frac{1}{k}\sum_{i=1}^k q_{(i)}\)
pvalues=rep(NULL, 778)
lfc=rep(NULL, 778)
bacteria = as.matrix(bacteria)
for(i in 1:778){
#1. Specify model
model=summary(glm(bacteria[i,]~biomass, family="poisson"))
#2. Compute p-values (and LFC while we're at it)
pvalues[i]=model$coefficient[2,4]
lfc[i]=model$coefficient[2,1]
}
## Warning: glm.fit: fitted rates numerically 0 occurred ## Warning: glm.fit: fitted rates numerically 0 occurred ## Warning: glm.fit: fitted rates numerically 0 occurred ## Warning: glm.fit: fitted rates numerically 0 occurred ## Warning: glm.fit: fitted rates numerically 0 occurred ## Warning: glm.fit: fitted rates numerically 0 occurred ## Warning: glm.fit: fitted rates numerically 0 occurred ## Warning: glm.fit: fitted rates numerically 0 occurred ## Warning: glm.fit: fitted rates numerically 0 occurred
Above gets p-values and log fold changes (more later) for log linear model
Yes warnings can and will happen in even simple generalized linear model
# 3. Estimate FDR and return TRUE if < .05 discoveries = p.adjust(pvalues, method = "BH")<.05 sum(discoveries)
## [1] 34
plot(pvalues, p.adjust(pvalues, method="BH"), ylab="FDR", xlim = c(0,.05),ylim=c(0,.5), col= as.numeric(discoveries)+1)
abline(h = 0.05,v=.05, lwd =c(5,1), col = c(2,1))
legend("topleft", inset=.05, col=2, pch=1, legend = "Discovery")
sum(pvalues<.05)
## [1] 135
778*0.05
## [1] 38.9
38.9/135
## [1] 0.2881481
Bonferroni procedure uses \(P<.05 / 778\) and controls Family-Wise Error Rate (FWER)
Q: How many discoveries would be made using Bonferroni rule?
sum(pvalues<.05/778)
## [1] 17
#1. Priors pi, beta estimated in Habiger, Watts, Anderson (2017)
pi = c(.69,.16,.15)
beta=c(0,-1.13,.78)
#2 Compute posterior:
p.0=exp(beta[1]*biomass)/sum(exp(beta[1]*biomass))
p.l=exp(beta[2]*biomass)/sum(exp(beta[2]*biomass))
p.u=exp(beta[3]*biomass)/sum(exp(beta[3]*biomass))
q<-function(y){ pi[1]*dmultinom(y,prob=p.0)/
(pi[1]*dmultinom(y,prob=p.0) + pi[2]*dmultinom(y,prob=p.l)+pi[3]*dmultinom(y, prob=p.u))}
qs=rep(NULL, 778)
for(i in 1:778){qs[i]=q(bacteria[i,])}
Notable differences in steps 1 and 2
Specified pi and beta values
Had to use Bayes theorem and my own code to compute \(q\) - not using glm()
#3 EBdiscoveries are ones with FDR<.05 s.q=sort(qs) FDR=cumsum(s.q)/1:778 k=sum(FDR<=.05) EBdiscoveries=qs<=s.q[k] sum(EBdiscoveries)
## [1] 97
plot(s.q, FDR, xlab="Post. Null", ylab = "FDR", xlim = c(0,.3), ylim = c(0, .1), col=I(FDR<.05)+1)
abline(h = 0.05, lwd =5, col = 2)
legend("topleft", inset=.05, col=2, pch=1, legend = "Discovery")
biplot(princomp(l.bacteria[,]),scale=0)
biplot(princomp(l.bacteria[discoveries,]),xlim = c(-3,8),ylim = c(-4, 4),scale=0)
biplot(princomp(l.bacteria[EBdiscoveries,]),xlim = c(-3,8),ylim = c(-4, 4), scale=0)
plot(lfc,p.adjust(pvalues, method="BH"), xlim = c(-3,3), xlab = "log fold change per gram", ylab=expression(hat(FDR)))
abline(h=0.05, lwd =3)
#index=sort(qs, index=TRUE)[[2]]
#points(lfc[index], FDR,pch=2, col=2)
legend("topright",legend=c("Classical","Emp Bayes"), pch = c(1,2), col = c(1,2))
plot(lfc,p.adjust(pvalues, method="BH"), xlim = c(-3,3), xlab = "log fold change per gram", ylab=expression(hat(FDR)))
abline(h=0.05, lwd =3)
index=sort(qs, index=TRUE)[[2]]
points(lfc[index], FDR,pch=2, col=2)
legend("topright",legend=c("Classical","Emp Bayes"), pch = c(1,2), col = c(1,2))
Properties: BH procedure controls the FDR If P-values are valid and, say, positively correlated.
Step 1. glm() function specified Poisson model and \(H_0:\) LFC\(_i=0\) (the log fold change per biomass gram for bacteria \(i\))
What about models for…
Over-dispersion: Try adding “family = quasipoisson” to glm(), or using nb.glm(), chisq.test(), lm() on logged data. Do you get the same results?
Zero-inflation: Recall lots of zeros for some bacteria
Normalization: Try glm(bacteria[i,]~biomass + colSums(bacteria),..). Do you get the same results? Why not median, etc.
Why \(H_0:\) LFC\(_i = 0\) and why not \(H_0:\) LFC\(_i\approx 0\)?
Step 2. Get \(Z_i = \widehat{\text{LFC}}_i/SE_i\) and \(P_i=\Pr_0(|Z_i|\geq |z_i|)\)
P-values are valid if we use the correct distribution in \(\Pr_0\).
Is Z\(\sim\) N(0,1) used in glm “correct”?
Over/under fitting and bias vs. variance:
We can prove standard errors (and p-values) go to 0 as total abundance increases…EVEN IF LFC\(_i = 0.000001\) in reality!
Do \(H_0:\) LFC\(_i\approx 0\) and \(H_0:\) LFC\(_i=0\) agree here?
Who cares?
Statistics Show Fred is Associated with Productivity!
Step 3. Compute \(\widehat{FDR}(p) = \frac{M p}{(\# p_i\leq p)}\) for each \(p\).
The numerator estimates # False Discoveries as M \(p\).
Can we use \(SE_i\) in weighted p-value procedures?
How can we choose weights?
How should criteria, data, etc. inform other weights
Properties: Empirical Bayes procedures should (asymptotically) control the FDR and make more discoveries than classical approach counterparts
Step 1. Assume 3 possible multinomial distributions and \(H_0:\)LFC\(_i=0\)
Parameter Estimation: Maximum likelihood estimators below. Other types?
Statistical Model: Why 3 multinomials? Why multinomial? These choices affect
Is \(H_0:\) LFC\(_i = 0\) or \(H_0:\) LFC\(_i\approx 0\) tested in above considerations?
Step 2. Estimate posterior null probability \(q_i = \Pr(LFC_i=0|\text{Bacteria}_i)\)
Are posterior null probabilities approximately “valid”?
Over/under fitting, over-dispersion, zero-inflation
Number of multinomials \(k\) “tunes” bias vs variance / over-under fitting
Is \(q_i\) the same for \(H_0:\) LFC\(_i\approx 0\) vs \(H_0:\) LFC\(_i=0\)?
Does number of multinomials (k=3, 5, 10,…) matter?
Does “\(\approx\)” matter?
Can Fred become “significant” for some \(k\) and “\(\approx\)” definitions?
Step 3. Estimate \(FDR(q_{(r)}) = \bar{q}\)
We could prove that ANY DISCOVERY RULE with \(\bar{q}\leq \alpha\) controls FDR!
Standard procedure maximize the number of discoveries
Quantity vs quality?
Why is the Benjamini and Hochber (1995) procedure one of the most highly cited statistics papers of all time?
Better than Bonferroni, especially when num tests is large
We can get usually get a p-value for quite complex models
Presume p-value is valid
When to use an Empirical Bayes methods?
You know how to select a good model and null hypotheses
You know how to estimate parameters well
| Deseq2, EdgeR | PMLfdr | |
|---|---|---|
| Approach | P-value | Empirical Bayes |
| Overdispersion | YES | YES + formula! |
| Over/Under Fitting | YES | YES + model selection |
| FDR Control/Efficiency | simul/no | proof+simul / proof+simul |
| Is Fred Productive | YES | NO |
| Misuse Potential | Low | Moderate |
With Prof Rudra and Mehedi Hassan (PhD 2028): Empirical Bayes Methods for false discovery rate control
Poisson mixtures vs. Multinomial mixtures?
Affect of normalization methods in above?
Include other covariates in model?
Victor Abadom (MS 2026):
When will \(\text{log}(Y+1) \sim N(\mu,\sigma^2)\) work well vs. \(Y\sim\) Negative Binomial, quasi-Poisson, etc.?
What about zero inflation?
Rhizosphere-Mediated Effects on Winter Wheat Yield in Oklahoma and the Development of Diagnostic Probes Targeting Beneficial Microorganisms. Ramos-Lopez, Espindola, et. al. (2026+)
Model: log_Abun\(_i\sim\beta_{0i} +\beta_{1i} f(\text{Phenotype}_{1i},..., \text{Phenotype}_{5i})\)
\(f()\) indicates high yield
12 trt groups (Soil x Variety x Feek Stage) + 3 sampling depths
HD replicability analysis could be used to identify bacteria differentially expressed in \(2\) + groups, \(3\) + groups, . . . ?
HD mediation analysis could be used to identify which bacteria are mediating relationships between factors and yield?
Challenges: EVERYTHING WE JUST TALKED ABOUT and
\[\text{logit}(\Pr(\text{Pathogen})) = \beta_1 EP_1 + \beta_2 EP_2 + ... +\beta_{1000} EP_{1000}\]
Sensitivity Analysis for FDR Estimation Under Misspecified Selection Model. Cao, Yi, Habiger (2026+)
Motivation
Working Selective Inference model for science
Question: Is selection mechanism really that simple? What if we’re wrong?