libraries

library("ggplot2")
library("lme4")
library("dplyr")
library("merTools")

Simulate data with fixed mean and variance

reps <- 50
samps <- 30
g_mean <- 0
g_var <- 4
Camden_OA <- NULL
for(i in seq_len(reps)){
  group_price <- rnorm(samps,g_mean,g_var)
  group_i <- data.frame(Trans_Log = group_price, MSOA11CD = paste0("group",i))
  Camden_OA <- rbind(Camden_OA, group_i)
}

Boxplot of groups

ggplot(Camden_OA, aes(x = MSOA11CD, y = Trans_Log)) +
  geom_boxplot() +
  theme_bw()

model fixed mean/variance groups

model.0 <- lmer(Trans_Log ~ 1 + (1|MSOA11CD), data = Camden_OA) 
summary(model.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Trans_Log ~ 1 + (1 | MSOA11CD)
##    Data: Camden_OA
## 
## REML criterion at convergence: 8405.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2469 -0.6442  0.0010  0.6773  2.7987 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  MSOA11CD (Intercept) 6.963e-15 8.344e-08
##  Residual             1.587e+01 3.984e+00
## Number of obs: 1500, groups:  MSOA11CD, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.03463    0.10286   0.337

diffence between intercept coeff and group mean

intercept_coeff <- coef(model.0)
group_mean <- group_by(Camden_OA, MSOA11CD) %>%
  summarise(mean = mean(Trans_Log))
plot_coef <- data.frame(coeff = intercept_coeff$MSOA11CD[,1], mean2 = group_mean$mean)
plot_coef$error <- abs((plot_coef$coeff - plot_coef$mean2)^2)
print(sum(plot_coef$error))
## [1] 22.30617
plot(plot_coef$coeff, plot_coef$mean2)

### Caterpillar plot

# Extract the random effect estimates:
MSOA_re <- REsim(model.0)
# Plot as a catipillar graph:
model.0_re <- plotREsim(MSOA_re)
# Check the results:
model.0_re

rerun same analysis with variable mean/variance

reps <- 50
samps <- 30
Camden_OA <- NULL
for(i in seq_len(reps)){
  g_mean <- sample(seq(0,5,0.1),1)
  g_var  <- sample(seq(0,2,0.1),1)
  group_price <- rnorm(samps,g_mean,g_var)
  group_i <- data.frame(Trans_Log = group_price, MSOA11CD = paste0("group",i))
  Camden_OA <- rbind(Camden_OA, group_i)
}

model.0 <- lmer(Trans_Log ~ 1 + (1|MSOA11CD), data = Camden_OA) 

intercept_coeff <- coef(model.0)
group_mean <- group_by(Camden_OA, MSOA11CD) %>%
  summarise(mean = mean(Trans_Log))
plot_coef <- data.frame(coeff = intercept_coeff$MSOA11CD[,1], mean2 = group_mean$mean)
plot_coef$error <- abs((plot_coef$coeff - plot_coef$mean2)^2)
sum(plot_coef$error)
## [1] 0.09691664
# ggplot(plot_coef, aes(x = coeff, y = mean2)) +
#   geom_line()

# Extract the random effect estimates:
MSOA_re <- REsim(model.0)
# Plot as a catipillar graph:
model.0_re <- plotREsim(MSOA_re)
# Check the results:
model.0_re