How to estimate variance for identical & fraternal twins with a mixed model

Question

I am using the twins data from OpenMx and am curious to estimate the empirical covariance structure for BMI differences in monozygous and dizygous twins:

library(OpenMx)
data(twinData)
twinData$dz <- as.numeric(twinData$zyg >= 3)

## empirical variance
by(twinData[, c('bmi1', 'bmi2')], twinData$dz, cov, use='complete.obs')
# twinData$dz: 0
# bmi1 bmi2
# bmi1 0.70 0.56
# bmi2 0.56 0.74
# ----------------------------------------------------------------- 
#   twinData$dz: 1
# bmi1 bmi2
# bmi1 0.91 0.41
# bmi2 0.41 0.88

To do this, I restructured the data to a long format with an indicator of dizygous twin group and twin pair ID.

twinData2 <- data.frame(bmi = with(twinData, c(bmi1, bmi2)))
twinData2$dz <- twinData$dz
twinData2$id <- factor(1:nrow(twinData))

Then I fit the following model and found:

summary(fit2 <- lmer(bmi ~ dz + (1 | id) + (0  + dz | id), data=twinData2)) 

...

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.4233   0.651   
 id.1     dz          0.0128   0.113   
 Residual             0.4207   0.649   
Number of obs: 7362, groups:  id, 3793

This predicts the BMI variance in dizygous twins to be 0.42 + 0.01 + 0.42 = 0.85 and in monozygous to be 0.84.

Is this calculation correct? Why are these predictions so different from the empirical variance?