# R/covjmcm_mcd.R In varjmcm: Estimations for the Covariance of Estimated Parameters in Joint Mean-Covariance Models

#### Documented in covjmcm_mcd

```#' Calculate the estimation of the covariance of estimated parameters in a MCD model, via the explicit formula.
#'
#' \code{covjmcm_mcd} gives an estimation of the covariance of estimated parameters in a MCD model using
#' the explicit formula, which is the inverse of the estimated Fisher's information matrix.
#'
#' @param object a fitted joint mean-covariance model of class "jmcmMod", returned by the function \code{jmcm}.
#' @return an estimated covariance matrix of the estimated parameters in a MCD model.
#' @examples
#' cattleA <- cattle[cattle\$group=='A', ]
#' fit.mcd <- jmcm(weight|id|I(ceiling(day/14+1))~1|1,
#'                data = cattleA, cov.method = "mcd",
#'                triple = c(8,3,4))
#' cov.mcd <- covjmcm_mcd(fit.mcd)
#' @references [1] Pourahmadi, M., "Maximum likelihood estimation of generalised linear models for multivariate normal covariance matrix," Biometrika
#' 87(2), 425–435 (2000).
#' @importFrom MASS ginv
#' @export
covjmcm_mcd <- function(object){
##object is a fitted jmcm model
if (missing(object)) stop("missing object.")
if(object@call\$cov.method!="mcd") stop("Method must be mcd")
m <- getJMCM(object, name="m")
n <- length(m)
q <- length(getJMCM(object,"gamma"))
p <- length(getJMCM(object,"beta"))
d <- length(getJMCM(object,"lambda"))
p1 <- p+d
p2 <- p+q+d

I <- matrix(0, p2, p2)

for(i in 1:n){

Xi <- getJMCM(object, name="X",sub.num=i)
Zi <- getJMCM(object, name="Z",sub.num=i)
Wi <- getJMCM(object, name="W",sub.num=i)
Sigma <- getJMCM(object, name="Sigma", sub.num=i)
T <- getJMCM(object, name="T", sub.num=i)
D <- as.vector(diag(getJMCM(object, name="D", sub.num=i)))

I11 <- t(Xi) %*% ginv(Sigma) %*% Xi

I22 <- (1/2)*(t(Zi)%*%Zi)

if(m[i]==1){I33 <- matrix(0,q,q)
I23 <- matrix(0,d,q)}

if(m[i]>1){
wi <- array(0, dim=c(q , q,  m[i]))
wi[,,1] <- matrix(0, q, q)
for(j in 2:m[i]){
for(k in 1:(j-1)){
for(l in 1:(j-1))
wi[,,j] <- wi[,,j] + (Sigma[k,l]/D[j])*
(tcrossprod(Wi[sum(1:(j-1))-(j-1)+k, ],Wi[sum(1:(j-1))-(j-1)+l, ]))
}
}

I33 <- apply(wi,c(1,2),sum)

A <- Sigma %*% t(T)
B <- matrix(0, m[i], q)

for (j in 2:m[i]) {
for (k in 1:(j-1))
B[j,] <- B[j,] + A[k,j] * Wi[sum(1:(j-1))-(j-1)+k, ]
}

I23 <- t(Zi) %*% ((1/D)*B)
}

I[1:p,1:p] <- I[1:p,1:p] + I11
I[(p+1):p1,(p+1):p1] <- I[(p+1):p1,(p+1):p1] + I22
I[(p1+1):p2,(p1+1):p2] <- I[(p1+1):p2,(p1+1):p2] + I33
I[(p+1):p1,(p1+1):p2] <- I[(p+1):p1,(p1+1):p2] + I23
}

I[(p1+1):p2,(p+1):p1] <- t(I[(p+1):p1,(p1+1):p2])
cov <- ginv(I)
cov
}
```

## Try the varjmcm package in your browser

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

varjmcm documentation built on April 14, 2020, 6:58 p.m.