Basicamente não existe segredo para reproduzir a mesma análise. O modo mais direto é declarar o modelo para estimar os parâmetros separados para cada nível do fator, depois montar as matrizes correspondentes às funções lineares e submetê-las à glht() ou qualquer outra função que permita inferências para funções lineares de parâmetros. O pacote car, contrast e gmodels têm funções equivalentes à glht(), minha preferida.

##-----------------------------------------------------------------------------

m0 <- lm(Y~0+estudo/X, DATA)
summary(m0)

m0$assign
tapply(coef(m0), m0$assign, mean)

## Matriz de pesos.
m <- rbind(rep(1, nlevels(DATA$estudo))/nlevels(DATA$estudo))

require(multcomp)

## Intercepto médio.
X <- cbind(m, 0*m)
summary(glht(m0, linfct=X), test=adjusted(type="none"))

## Inclinação média.
X <- cbind(0*m, m)
summary(glht(m0, linfct=X), test=adjusted(type="none"))

## As inclinações para cada nível saem no próprio summary() do modelo
## quando se declara o modelo na forma ~0+fator/numérica. Para
## conhecimento, pode ser feito assim.

M <- diag(nlevels(DATA$estudo))
X <- cbind(M, 0*M)
summary(glht(m0, linfct=X), test=adjusted(type="fdr"))

X <- cbind(0*M, M)
summary(glht(m0, linfct=X), test=adjusted(type="fdr"))

##-----------------------------------------------------------------------------

À disposição.

Walmes.