library(nlme)

# Data prep
dat <- clubSandwich::MortalityRates[, c("year", "state", "mrate", "cause")]
dat <- dat[dat$cause == "Motor Vehicle" & dat$year %in% 1970:1975, ]
dat <- dat[, c("year", "state", "mrate")]
str(dat)
'data.frame':   306 obs. of  3 variables:
 $ year : int  1970 1971 1972 1973 1974 1975 1970 1971 1972 1973 ...
 $ state: int  1 1 1 1 1 1 2 2 2 2 ...
 $ mrate: num  62.2 62.4 73.6 63.4 60.3 ...

dat$year.f <- factor(dat$year) # Make time dummies
dat$state.f <- factor(dat$state) # Make state dummies

(m0 <- lme(mrate ~ year.f, dat, ~ 1 | state.f))
Linear mixed-effects model fit by REML
  Data: dat
  Log-restricted-likelihood: -1231.598
  Fixed: mrate ~ year.f
(Intercept)  year.f1971  year.f1972  year.f1973  year.f1974  year.f1975
 66.3498221  -3.4835741   0.2247485   0.3086386  -8.1646117 -10.6781779

Random effects:
 Formula: ~1 | state.f
        (Intercept) Residual
StdDev:    22.29717  10.7256

Number of Observations: 306
Number of Groups: 51

as.numeric(VarCorr(m0))[1] / (
  as.numeric(VarCorr(m0))[1] + var(predict(m0, level = 0)) + sigma(m0) ** 2)
[1] 0.7876077

mg0 <- gls(mrate ~ year.f, dat, corCompSymm(form = ~ 1 | state.f))
mg1 <- gls(mrate ~ year.f, dat) # A regular regression model ignoring states
anova(mg0, mg1) # Model comparison

    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
mg0     1  8 2479.197 2508.827 -1231.598
mg1     2  7 2814.073 2839.999 -1400.036 1 vs 2 336.8762  <.0001