```

References

set.seed(102)

## 4-category response variable is defined
probs <- c(0.1,0.3,0.5,0.2)
X <- rmultinom(n = 100, size = 1, prob = probs)
dat <- data.frame(X  =  as.vector(seq_along(probs) %*% X),
                  ## Group variable produced independently
                  group = c(rep(1,20), rep(2,80)))
## But groups are nonetheless somewhat imbalanced
datP <- prop.table(xtabs( ~ X + group, dat), margin = 2)

## Treated group vector (dropping the first outcome category)
(T <- datP[-1,1])
##   2   3   4 
## 0.3 0.3 0.3
## Control group (dropping the first outcome category)
(C <- datP[-1,2])
##      2      3      4 
## 0.2750 0.4875 0.1625
## Mean difference vector (T - C)
(meanDiffVector <- (T - C))
##       2       3       4 
##  0.0250 -0.1875  0.1375
## Variance-covariance matrix has to be defined for both groups
## As defined in reference (2) and used in reference (1)
## Non-diagonal elements are negated cross product of k-th element and l-th element
vcovT <- -1 * outer(T, T)
## Diagonal elements where k = l are p(1-p) similarly to binary case
diag(vcovT) <- T * (1-T)
vcovT
##       2     3     4
## 2  0.21 -0.09 -0.09
## 3 -0.09  0.21 -0.09
## 4 -0.09 -0.09  0.21
## similarly for C vector
vcovC <- -1 * outer(C, C)
diag(vcovC) <- C * (1-C)
vcovC
##            2           3           4
## 2  0.1993750 -0.13406250 -0.04468750
## 3 -0.1340625  0.24984375 -0.07921875
## 4 -0.0446875 -0.07921875  0.13609375
## These variance covariance matrices together constitute S matrix in reference (1)
## Division by 2 is element-wise division
(S <- (vcovT + vcovC) / 2)
##             2           3           4
## 2  0.20468750 -0.11203125 -0.06734375
## 3 -0.11203125  0.22992187 -0.08460938
## 4 -0.06734375 -0.08460938  0.17304687
## Invert S using solve() or MASS::ginv() (generalized inverse)
(Sinv <- solve(S))
##          2        3        4
## 2 14.91717 11.46822 11.41250
## 3 11.46822 14.12023 11.36696
## 4 11.41250 11.36696 15.77787
## (T - C)^T S^(-1) (T - C)
(smdCat <- drop(t(meanDiffVector) %*% Sinv %*% t(t(meanDiffVector))))
## [1] 0.1888756

The same can be done using the cov() function a dummy variable matrix, but the (N - 1) denominator has to be adjusted to be the same.

## Append dummy variables
dat <- cbind(dat, dummies::dummy(dat$X))
names(dat)[3:6] <- paste0("X",1:4)

## Regular vcov using cov() on 3 of dummies
(vcovT2 <- 19/20 * cov(subset(dat, group == 1)[c("X2","X3","X4")]))
##       X2    X3    X4
## X2  0.21 -0.09 -0.09
## X3 -0.09  0.21 -0.09
## X4 -0.09 -0.09  0.21
(vcovC2 <- 79/80 * cov(subset(dat, group == 2)[c("X2","X3","X4")]))
##            X2          X3          X4
## X2  0.1993750 -0.13406250 -0.04468750
## X3 -0.1340625  0.24984375 -0.07921875
## X4 -0.0446875 -0.07921875  0.13609375
## The rest is the same
(S2 <- (vcovT2 + vcovC2) / 2)
##             X2          X3          X4
## X2  0.20468750 -0.11203125 -0.06734375
## X3 -0.11203125  0.22992187 -0.08460938
## X4 -0.06734375 -0.08460938  0.17304687
(Sinv2 <- solve(S2))
##          X2       X3       X4
## X2 14.91717 11.46822 11.41250
## X3 11.46822 14.12023 11.36696
## X4 11.41250 11.36696 15.77787
(smdCat2 <- drop(t(meanDiffVector) %*% Sinv2 %*% t(t(meanDiffVector))))
## [1] 0.1888756

Comparison

list(vcovT = vcovT, vcovT2 = vcovT2)
## $vcovT
##       2     3     4
## 2  0.21 -0.09 -0.09
## 3 -0.09  0.21 -0.09
## 4 -0.09 -0.09  0.21
## 
## $vcovT2
##       X2    X3    X4
## X2  0.21 -0.09 -0.09
## X3 -0.09  0.21 -0.09
## X4 -0.09 -0.09  0.21
list(vcovC = vcovC, vcovC2 = vcovC2)
## $vcovC
##            2           3           4
## 2  0.1993750 -0.13406250 -0.04468750
## 3 -0.1340625  0.24984375 -0.07921875
## 4 -0.0446875 -0.07921875  0.13609375
## 
## $vcovC2
##            X2          X3          X4
## X2  0.1993750 -0.13406250 -0.04468750
## X3 -0.1340625  0.24984375 -0.07921875
## X4 -0.0446875 -0.07921875  0.13609375
list(S = S, S2 = S2)
## $S
##             2           3           4
## 2  0.20468750 -0.11203125 -0.06734375
## 3 -0.11203125  0.22992187 -0.08460938
## 4 -0.06734375 -0.08460938  0.17304687
## 
## $S2
##             X2          X3          X4
## X2  0.20468750 -0.11203125 -0.06734375
## X3 -0.11203125  0.22992187 -0.08460938
## X4 -0.06734375 -0.08460938  0.17304687
list(Sinv = Sinv, Sinv2 = Sinv2)
## $Sinv
##          2        3        4
## 2 14.91717 11.46822 11.41250
## 3 11.46822 14.12023 11.36696
## 4 11.41250 11.36696 15.77787
## 
## $Sinv2
##          X2       X3       X4
## X2 14.91717 11.46822 11.41250
## X3 11.46822 14.12023 11.36696
## X4 11.41250 11.36696 15.77787
list(smdCat = smdCat, smdCat2 = smdCat2)
## $smdCat
## [1] 0.1888756
## 
## $smdCat2
## [1] 0.1888756