controlling for clustering at id level in mixed effects model I have one group ($n=40$) of subjects pre- and post-tested (time; coded $0$ and $1$) on a continuous variable (y). I also have a gender variable coded $0$ and $1$.
Using R, I was wondering how I could run a mixed effects model using time point and gender as the independent variables, controlling for clustering at id level?
P.S. I'm interested in detecting a meaningful difference among the average scores between the pre- and post-scores as a whole.
Here is my data in R:
dat <- {
 set.seed(62)
 data.frame(id = rep(1:40, each = 2), y = rnorm(80, 20, 5), gender = rep(rbinom(40, 1, .5), each = 2), 
 time = rep(0:1, 40))

}

 A: This is a pretty standard specification. 


*

*Assuming that the conditional distribution of y (i.e., residuals from the model) are reasonably close to Gaussian, a linear mixed model is appropriate.

*Interacting fixed effects of gender and time and a random intercept varying across id should do what you want. 


*

*It doesn't make sense to test for differences in the effect gender by id because [unless you're working with people who are gender transitioning, or certain species of sex-transitioning fish] each individual will only have a single id. 

*The variation in time effect by id is subsumed in the residual variation because there are only two time observations per person.


*note that the interpretation of the main effects of gender and time will depend on the contrasts used; switch from the default treatment contrasts to sum-to-zero contrasts if you want average effects. (Or use the emmeans package to do this post hoc.)


Frequentist
## oldest, provides p-values by default
nlme::lme(y~gender*time, random= ~1|id, data=dat)
## faster, more flexible; load the lmerTest package
## instead of lme4 to get p-values
lme4::lmer(y~gender*time+(1|id),data=dat)
## recent, even more flexible
glmmTMB::glmmTMB(y~gender*time+(1|id),data=dat)

Bayesian
Since you used the [hierarchical-bayesian] tag:
MCMCglmm::MCMCglmm(y~gender*time, random=~id, data=dat)
brms::brm(y~gender*time+(1|id),data=dat)

A: Ben has provided a great answer.  Another approach would be to use the gee or geepack packages for marginal model (using the default gaussian family for your data).  You would have to pick an appropriate covariance structure and specify that in the corstr argument.  In your case this would be: 
geepack::geeglm(y~time*gender, id=id, data=dat, corstr="ar1", waves=time)

or 
gee::gee(formula = y ~ time * gender, id = id, data = dat, corstr = "unstructured")

