R code to fit varying intercepts with nested structure linear model

## sim some data and fit an varying intercept model
##

a <- 1:3
b <- 1:4

intercepts <- rnorm(12) 
intnames <- outer(a, b, FUN = function(x, y) { paste(x, y, sep = ':')} )
names(intercepts) <- as.vector(t(intnames)) 

N <- 10000; 
x <- round(abs(rnorm(N)) * 10, 2)
noise <- rnorm(N); 
ai <- sample(a, N, replace = TRUE) 
bi <- sample(b, N, replace = TRUE) 
beta <- 4; 
y <- intercepts[ai * 4 + bi - 4] + x * beta + noise 

ai <- as.factor(ai)
bi <- as.factor(bi)
table(ai:bi) 
print(intercepts) 
# varying intercepts model specification
fita <- lm(y ~ -1 + ai : bi + x) 
# in R lm and lmer's model specification, A * B = A + B + A : B
# where A : B is the interaction of A and B. So in the above,
# ai : bi is used to specify varying intercepts and -1 to
# get rid of the global intercept.