mixed: Estimating parameters of a linear mixed model (LMM) with a...

In gajdosandrej/fdslrm: R package for modeling and prediction of time series using linear mixed models

Description

mixed(y, X, Z, dim, s20, method, lambda, adaptRW) computes ML, REML, MINQE(I), MINQE(U,I), BLUE(b), BLUP(u) by Henderson's Mixed Model Equations Algorithm. This is R version of the original Matlab program created by Viktor Witkovsky (Witkovsky, 2000).

Model: Y=X*b+Z*u+e,

b=(b_1',...,b_f')' and u=(u_1',...,u_r')', E(u)=0, Var(u)=diag(sigma^2_i*I_{m_i}), i=1,...,r E(e)=0, Var(e)=sigma^2_{r+1}*I_n, Var(y)=Sig=sum_{i=1}^{r+1} sigma^2_i*Sig_i. We assume normality and independence of u and e.

Usage

mixed(y, X, Z, dim, s20, method, lambda, adaptRW)

Arguments

y	n-dimensional vector of observations.
X	(n * k)-design matrix for fixed effects b=(b_1,...,b_f), typically X=[X_1,...,X_f] for some X_i.
Z	(n * m)-design matrix for random efects u=(u_1,...,u_r), typically Z=[Z_1,...,Z_r] for some Z_i.
dim	vector of dimensions of u_i, i=1,...,r, dim=(m_1,...,m_r), m=sum(dim).
s20	a prior choice of the variance components, s20=(s20_1,...,s20_r,s20_r+1). SHOULD BE POSITIVE for method>0.
method	method of estimation of variance components; 0:NO estimation, 1:ML, 2:REML, 3:MINQE(I), 4:MINQE(U,I)
lambda	regularization parameter used for ridge regression weights, default value is lambda = numeric() (use standatd estimation procedure, i.e. no regularized ridge estimation procedure).
adaptRW	flag for using adaptive method for the ridge matrix weights. adaptRW = FALSE, the used ridge matrix is Rw = lambda * diag(rep(1,k)). If adaptRW = TRUE, Rw = lambda * diag(weights), where weights = rep(1,k)/abs(b) and b is fixed effect estimate in current iteration. Default value is adaptRW = FALSE.

Value

List with the following components:

• s2 - estimated vector of variance components (sigma^2_1,..., sigma^2_{r+1})'. A warning message appears if some of the estimated variance components is negative or equal to zero. In such cases the calculated Fisher information matrices are inadequate.

• b - k-dimensional vector of estimated fixed effects beta, b=(b_1,...,b_f)=(X'Sig^{-1}X)^{+}X'Sig^{-1}y.

• u - m-dimensional vector of EBLUP's of random effects U, u=(u_1,...,u_r).

• Is2 - Fisher information matrix for variance components; if method=0 the output is empty matrix Is2; if metod=3 or method=4 the output is inversion of the covariance matrix of MINQE calculated at estimated s2.

• C - g-inverse of Henderson's MME matrix, where C=MPinv([XX XZ; XZ' ZZ+inv(D)*s0]/s0), if inv(D) exists or C=s0*[I 0; 0 D]*pinv([XX XZ*D; XZ' V]) otherwise

• H - Criterial matrix for MINQE calculated at priors s20; if method=3: H_ij=trace(Sig_0^{-1}*Sig_i*Sig_0^{-1}*Sig_j), if method=4: H_ij=trace((M*Sig_0*M)^{+}*Sig_i*(M*Sig_0*M)^{+}*Sig_j)

• q - (r+1)-dimensional vector of MINQE(U,I) quadratic forms calculated at prior values s20; if method=0,1,2 the output is empty vector q, otherwise q_i=y'*(M*Sig_0*M)^{+}*Sig_i*(M*Sig_0*M)^{+}*y.

• loglik - Log-likelihood evaluated at the estimated parameters; if method=1: loglik=log-likelihood(ML), if method=2: loglik=log-likelihood(REML), if method=3: or method=4 loglik=numeric(), if method=0: loglik=informative value of log of the joint pdf of (y,u).

• loops - Number of loops.

Note

Ver.: 23-Apr-2020 19:44:40.

References

Witkovsky, V.: MATLAB Algorithm for solving Henderson's Mixed Model Equations. Technical Report No. 3/2000, Institute of Measurement Science, Slovak Academy of Sciences, Bratislava, 2000. https://www.mathworks.com/matlabcentral/fileexchange/200-mixed?s_tid=prof_contriblnk

Searle, S.R., Cassela, G., McCulloch, C.E.: Variance Components. John Wiley & Sons, INC., New York, 1992. (pp. 275-286).

Examples

## EXAMPLE 1:
IncN <- t(matrix(c(4, 5, 8, 9, 5, 10, 15, 20), 4, 2))
matrices <- Design2(IncN)
n <- nrow(matrices$A)
n1 <- ncol(matrices$B)
n2 <- ncol(matrices$C)
X <- cbind(matrix(1, n, 1), matrices$A)
Z <- cbind(matrices$B, matrices$C)
btrue <- Conj(c(1, 2, 3))
s2true <- c(0.5, 3, 1)
u1 <- sqrt(s2true[1]) * rnorm(n1)
u2 <- sqrt(s2true[2]) * rnorm(n2)
u <- c(u1, u2)
e <- sqrt(s2true[3]) * rnorm(n)
y <- as.vector(X %*% btrue + Z %*% u + e)
dim <- c(n1, n2)
s20 <- c(1, 1, 1)
method <- 2     # 0:NONE, 1:ML, 2:REML, 3:MINQE(I), 4:MINQE(U,I)
result1 <- mixed(y, X, Z, dim, s20, method)
result1$s2
result1$b
result1$u
result1$Is2
result1$C
result1$H
result1$q
result1$loglik
result1$loops

gajdosandrej/fdslrm documentation built on April 28, 2020, 11:35 a.m.