Estimating and comparing treatment effects with participants that have repeated measurments Suppose we run an experiment with 50 participants. We take repeated measurements from each participant: each undergoes three different treatments during which we measure their reactions (our continuous target variable). A toy example of such dataset could be generated as follows:
# Example Data
treat_1 = rnorm(50,10,2)
treat_2 = rnorm(50,20,2)
treat_3 = rnorm(50,20,2)

df = data.frame(participant_id = rep(1:50,3), 
                treat_type = c(rep("1",50), rep("2",50), rep("3",50)),
                y = c(treat_1,treat_2,treat_3))

I am seeing people saying that when repeated measurements are involved we cannot simply run a linear regression and we need more sophisticated approaches such as mixed models and etc. But personally I don't understand what can go wrong with the following approach:
fit = lm(y ~ -1 + treat_type, data = df)

The only issue I can see, is that we are not accounting for heteroscedasticity (which is not the case in this example).
My question is if the approach I just showed is OK for estimating the treatment effects on the target variable considering that our dataset includes repeated measurements (three observation per each participant)? If not can someone please explain the reasons?
As a side question, why bother using t-test or ANOVA (which are common for comparing groups) that do not even provide any distributional information on the effects and just give p-values on whether there is a significant difference among the groups when we can fit a model that does give such information?
 A: That is a perfectly fine way to estimate the treatment effect. You can increase your precision by including a fixed effect for participant, i.e.,
fit = lm(y ~ -1 + treat_type + factor(participant_id), data = df)

This will not change the estimate if all participants have an equal number of observations for each treatment, but it will reduce the residual variance and increase precision (same reason we use a paired-samples t-test instead of an independent samples t-test when we have two treatment levels). This is equivalent to a repeated measures ANOVA.
You can also include participant ID as a random effect, e.g.,
fit = lme4::lmer(y ~ treat_type + (1|participant_id), data = df)

This should yield the same estimate as the fixed effect approach but the precision might be a little different.
Note that you should account for possible heteroscedasticity by using a robust standard error, which is easier with a fixed effects model. Some might note that observations from the same participant are correlated with each other and so you should use a cluster-robust standard error. This is not necessary as explained by Abadie et al. (2022). The usual robust standard error is sufficient.
I would recommend using estimatr::lm_robust() to do this, e.g.,
fit = estimatr::lm_robust(y ~ treat_type, data = df,
                          fixed_effects = ~participant_id)

This automatically incorporates the fixed effects and produces robust (non-clustered) standard errors.
