The package
A difficulty in calculating the moments of a multivariate truncated normal distribution is that the marginal distributions from a truncated normal distribution are not truncated normal. Numerous studies have been performed for the calculations under different conditions, including the type of truncation (one-sided or double-sided) and number of variables (univariate, bivariate, or multivariate). However, as of this writing, there are two R packages related to the moment calculation. The first is based on the work of Manjunath and Wilhelm (2012). They provided explicit representations for the mean and variance of a multivariate normal distribution with arbitrary double truncation using the moment generating function. These results were implemented in
One advantage of the
The general form of the model equation for a mixed model is,
where
The
Let
where
The evaluation of Pr[
The
The moment calculation for a truncated univariate normal distribution is implemented in the utrunmnt() function:
utrunmnt(kappa, mu = 0., sd = 1., lower = -Inf, upper = Inf)
Here, the kappa argument specifies the order of moment and must be provided as a nonnegative integer. The other arguments are clear by name. For example, if a normal distribution with mean 5 and standard deviation 1 is upper truncated at 10, the fourth moment is computed with the command:
> utrunmnt(4, mu = 5, sd = 1, upper = 10)
[1] 777.9971
The following R (R Core Team, 2019) code mimics the situation of a mixed-effects probit regression model. It generate a 5-dimensional latent variable from , where
set.seed(123)
sigma2e <- 1
sigma2a <- 2
n <- 5
mu <- seq(-1,1, length.out = n)
y <- mu + rnorm(1, sd = sqrt(sigma2a)) + rnorm(n, sd = sqrt(sigma2e))
y <- ifelse( y > 0, 1, 0)
a <- rep(-Inf, n)
b <- rep(Inf, n)
a[y >= 0] <- 0
b[y < 0] <- 0
To calculate the moment, you must first create an mtrunmnt object using mtrunmnt():
mtrunmnt(mu, Sigma2 = 1, lower = rep(-Inf, length(mu)), upper = rep(Inf, length(mu)),
Z = matrix(rep(1, length(mu)), ncol = 1), D = matrix(1, ncol = 1, nrow = 1),
nGH = 32)
The Sigma2 argument specifies the variances of the error terms and must be specified as a scalar, or a vector of length equal to mu. Given a scalar, the variances are assumed to be equal. The design matrix and variance-covariance matrix of the random effect are passed in the Z and D arguments, respectively. By default, the Z argument is provided as a vector of ones of the same length as mu, and should be changed accordingly when necessary. If the random components are independent and the variances are the same, you can pass the variance as a scalar to D. The nGH argument specifying the number of quadrature points is set to 32. The other arguments are self-explanatory.
The previously generated artificial data can be passed to the mtrunmnt() to create a mtrunmnt object. Once you have created the object, pass it to probntrun(), prodmnt() and meanvar() as follows:
> tobj <- mtrunmnt(mu, Sigma2 = sigma2e, lower = a, upper = b, D = sigma2a)
> probntrun(tobj)
[1] 0.01453587
> prodmnt(tobj, c(2,2,0,0,0))
[1] 1.612259
> c(prodmnt(tobj, c(4,0,0,0,0)), prodmnt(tobj, c(0,4,0,0,0)))
[1] 21.433824 1.482977
> > meanvar(tobj)
$tmean
[1] -1.4996487 0.5993896 -0.9354669 -0.7448227 1.1792822
$tvar
[,1] [,2] [,3] [,4] [,5]
[1,] 0.84907923 0.05062909 0.09792397 0.07408044 0.11730096
[2,] 0.05062909 0.26332061 0.03141038 0.02362801 0.04374713
[3,] 0.09792397 0.03141038 0.49954827 0.04749122 0.07292303
[4,] 0.07408044 0.02362801 0.04749122 0.36549384 0.05487201
[5,] 0.11730096 0.04374713 0.07292303 0.05487201 0.66073805
The probntrun() and meanvar() have only one argument specifying a trunmnt object, and compute Pr[
As we have seen, the
It is well known that the Gaussian-Hermite quadrature method gives an exact integration if integral function is a polynomial of order less than or equal to 2
Thus, it is easy to see that the functions enclosed in parentheses in
We want to verify the accuracy of
The Gauss-Hermit quadrature method is known to be efficient for low-dimensional integration. Therefore, the preceding example performed under the best conditions for trunmnt. In fact, it works best with random-effects one-way layout models because it’s essentially a one-dimensional integration problem because of the independence.
Consider an artificial example. Suppose
Z <- model.matrix( ~ factor(rep(c(0,1), each = 5)) - 1)
D <- diag(c(sigma2a, sigma2a))
mu2 <- c(mu, mu)
a2 <- c(a,a)
b2 <- c(b,b)
tobj2 <- mtrunmnt(mu2, Sigma2=sigma2e, lower = a2, upper = b2, Z = Z, D = D)
meanvar(tobj2)
We do not give the means and variance-covariance matrices of
$tvar
[,1] [,2] [,3] [,4] [,5]
[6,] 2.664535e-15 -3.330669e-16 8.881784e-16 1.110223e-15 -2.220446e-15
[7,] -1.110223e-15 4.996004e-16 -3.330669e-16 -3.885781e-16 5.551115e-16
[8,] 1.998401e-15 -4.440892e-16 1.110223e-16 3.330669e-16 -8.881784e-16
[9,] 1.998401e-15 -2.775558e-16 8.881784e-16 4.440892e-16 -8.881784e-16
[10,] -1.776357e-15 2.220446e-16 -6.661338e-16 -2.220446e-16 8.881784e-16
It’s worth noting that the
This paper introduces the
Summary of the functions in the
Function | Description |
meanvar | Calculate the mean and variance |
mtrunmnt | Create a mtrunmnt object with given parameters |
prodmnt | Calculate the product moment |
probntrun | Calculate the probability of not being truncated |
utrunmnt | Calculate the moment of univariate truncated normal distribution |