# R/Ddf.R In MissMech: Testing Homoscedasticity, Multivariate Normality, and Missing Completely at Random

```Ddf <- function(data, mu, sig) {
y <- data
n <- nrow(y)
p <- ncol(y)
ns <- p * (p + 1) / 2
nparam <- ns + p
ddss <- matrix(0, ns, ns)
ddmm <- matrix(0, p, p)
ddsm <- matrix(0, ns, p)
for (i in 1:n) {
obs <- which(!is.na(y[i, ]))
lo <- length(obs)
tmp <- cbind(rep(obs, 1, each = lo), rep(obs, lo))
tmp <- matrix(tmp[tmp[, 1] >= tmp[, 2], ], ncol = 2)
lolo = lo*(lo+1)/2
subsig <- sig[obs, obs]
submu <- mu[obs, ]
temp <- matrix(y[i, obs] - submu, nrow = 1)
a <- solve(subsig)
b <- a %*% (2 * t(temp) %*% temp - subsig) * a
d <- temp %*% a
ddimm <- 2 * a
#==================== DD(mu, mu)
ddmm[obs, obs] <- ddmm[obs, obs] + ddimm
#==================== DD(sig, mu)
rcnt <- 0
ddism <- matrix(0, lolo, lo)
for (k in 1:lo) {
for (l in 1:k) {
rcnt <- rcnt + 1
ccnt <- 0
for (kk in 1:lo) {
ccnt <- ccnt + 1
ddism[rcnt, ccnt] <- 2 * (1 - 0.5 * (k == l)) *
(a[kk, l] %*% d[k] + a[kk, k] %*% d[l])
}
}
}
for (k in 1:lolo) {
par1 <- tmp[k, 1] * (tmp[k, 1] -1) / 2 + tmp[k, 2]
for (j in 1:lo) {
ddsm[par1, obs[j]] <- ddsm[par1, obs[j]] + ddism[k, j]
}
}
#------------------------ Test part
ssi <- matrix(0, lolo, lolo)
for (m in 1:lolo) {
u <- which(obs == tmp[m, 1])
v <- which(obs == tmp[m, 2])
for (q in 1:m) {
k <- which(obs == tmp[q, 1])
l <- which(obs == tmp[q, 2])
ssi[m, q] <- (b[v, k] * a[l, u] + b[v, l] * a[k, u] + b[u, k] * a [l, v] +
b[u, l] * a[k, v]) * (1 - 0.5 * (u == v)) * (1 - 0.5 * (k == l))
}
}
for (k in 1:lolo) {
par1 <- tmp[k, 1] * (tmp[k, 1] -1) / 2 + tmp[k, 2]
for (l in 1:k) {
par2 <- tmp[l, 1] * (tmp[l, 1] - 1) / 2 + tmp[l, 2]
ddss[par1, par2] <- ddss[par1, par2] + ssi[k, l]
ddss[par2, par1] <- ddss[par1, par2]
}
}
}
dd <- -1 * rbind(cbind(ddmm, t(ddsm)), cbind(ddsm, ddss)) / 2
se <- -solve(dd)
list(dd = dd, se = se)
}
```

## Try the MissMech package in your browser

Any scripts or data that you put into this service are public.

MissMech documentation built on May 2, 2019, 1:08 p.m.