Luiz,

Você terá que mudar a função que está usando.

library(nlme)

table(x)

lattice::xyplot(y ~ x)

modelo0 <- aov(y ~ x)
a <- anova(modelo0)
a

# Componentes de variância pelo método ANOVA.
c("sigma^2_b" = (a["x", "Mean Sq"] - a["Residuals", "Mean Sq"])/6,
  "sigma^2_e" = a["Residuals", "Mean Sq"])

modelo1 <- lme(y ~ 1,
               random = ~1 | x,
               method = "REML")
modelo1

# Componentes de variância pelo método REML.
VarCorr(modelo1)

modelo2 <- lme(y ~ 1,
               random = ~1 | x,
               weights = varExp(),
               method = "ML")
modelo2

# Componentes de variância pelo método REML.
# ATTENTION: não são imediatamente comparáveis com os anteriores.
VarCorr(modelo2)

À disposição.
Walmes.