Group vs individual data: how to correct standard errors?

I have a set of 5 variables: Y is a response, and X1, X2, x3, and X4 are the predictors. There are 3,000 observations, and I fit a regression model in R:

fit <- lm(Y ~ X1 + X2 +X3 + x4, data = mydata)

The problem is that my data relate to 3,000 small units (parishes) nested in 40 larger units (counties), and Y, X1, and X2 are measured at the parish level, whilst X3 and X4 are only available at a county level. I am concerned that the lm function will underestimate the standard errors for the coefficients for the last two variables. Is there a way to correct this?

Mixed models! Something like

library(lmerTest)  ## so you can get p-values on fixed effects
fit <- lmer(Y ~ X1 + X2 + X3 + X4 + (1 + X1 + X2 | county), data = mydata)

The random-effect term allows the intercept, as well as the parish-level effect of X1 and X2, to vary across counties.

The logic here is that the fixed effects X1 + x2 + X3 + X4 specify that the model should allow these variables to affect the response at the population level (i.e., overall), while the random effects allow for the intercepts and the effects of X1 and X2 to vary across counties. The general rules are that (1) we can include any variable in a random-effects term that varies within the levels of the grouping variables (i.e., since X1 and X2 vary at the parish level, we can estimate the variance of their effects across counties; since X3 and X4 vary only between and not within counties, we cannot estimate the variation of their effects across counties) and (2) in general, because the random effects are zero-centered, we should usually include a fixed effect corresponding to each random-effects term.

For what it's worth you can use the equatiomatic package to extract LaTeX-formatted model specifications, see the vignette.