Hi @SingeunOh,

Absolutely you can use linear models or GLMs with categorical data. In fact, the common ANOVA is essentially a special case of a linear model. This is especially pretty straightforward to do with alpha diversity because your response variable would be numeric and often follows a normal distribution.

Below is a typical workflow I do myself in R with alpha diversity that may be helpful for you. It assumes you have a table that has rows as sites (ie. samples, participants, etc.) and meta data in your columns (ex. season, infection, richness, etc.)

```
#required libraries
library(tidyverse)
library(sjPlot)
library(car)
#build basic model
mod1 <- lm(alpha ~ season + infection, data=mydata)
#run a marginal test Anova (Type II)
mod1 %>% car::Anova() #optional, add test.statistic="F" for F-statistic
#you could also just use default aov() or summary() in R instead
#of the car package version, but that will be a sequential test by default
#check the model fit
par(mfrow=c(2,2))
sjPlot::plot_model(mod1, type="diag") %>%
sjPlot::plot_grid(.)
#if model diagnostics don't look good, you can try :
# - transforming the data first
# - Use glm() instead of lm() with a different link function (i.e. binomial, poisson, log, logit, etc.)
# - If nothing fits, then perhaps just stick with Wilcoxon and explain the above process
#get a nice table summary of the model results
sjPlot::tab_model(mod1)
#plot model results
sjPlot::plot_model(mod1, sort.est = TRUE, show.values = TRUE)
#if you want the exact values from this plot,
sjPlot::plot_model(mod1, sort.est = TRUE, show.values = TRUE) %>%
.$data %>%
tibble()
```

Also, within QIIME 2, you can use the q2-longitudinal plugin's `linear-mixed-effects`

action that will take R style formulas. Though this won't have flexibility to do any transformation or a glm, and won't get in depth model diagnostics.

As a note though, I would strongly recommend either reading a bit about linear models first or chat with a statistician about your specific project goals to make sure you are in fact addressing your specific questions, you're well powered, and the model you run is a decent fit. For ex. if you need to test the interaction term between season and an infection status you would change your formula to `alpha ~ season*infection`

which will demand different and more complex interpretation of results.

Hope that helps