R/BBmm_VarEst.R

In PROreg: Patient Reported Outcomes Regression Analysis

VarEst <- function(y,m,p,X,Z,u,nRand,nComp,nRandComp,OLDall.sigma,OLDphi,q,nDim){

  # Number of observations
  nObs <- length(y)/nDim

  # The wiegth matrix S
  s <- c(p*(1-p))

  # Defining the variance of the random effects
  d <- d. <- NULL
  for (i in 1:nComp){
    d <- c(d,rep(OLDall.sigma[i]^2,nRandComp[i]))
    d. <- c(d.,rep(1/(OLDall.sigma[i])^2,nRandComp[i]))
  }
  D <- diag(d)
  D. <- diag(d.)

  #-----------------------------------------------------------------------------
  # ESTIMATION OF PSI
  #-----------------------------------------------------------------------------
  function.psi <- function(psi){
    # Compute j and v
    v <- NULL
    hpsi <- NULL
    for (i in 1:nDim){
      for (l in ((i-1)*nObs+1):(i*nObs)){
        v1 <- 0
        v2 <- 0
        hpsi1 <- 0
        hpsi2 <- 0
        hpsi3 <- 0

        if (y[l]==0){
          v1 <- 0
          hpsi1 <- 0
        }else{
          for (k in 0:(y[l]-1)){
            v1 <- v1+1/(p[l]+k*exp(psi[i]))^2
            hpsi1 <- hpsi1+log(p[l]+k*exp(psi[i]))
          }
        }
        if (y[l]==m[l]){
          v2 <- 0
          hpsi2 <- 0
        }else{
          for (k in 0:(m[l]-y[l]-1)){
            v2 <- v2+1/(1-p[l]+k*exp(psi[i]))^2
            hpsi2 <- hpsi2+log(1-p[l]+k*exp(psi[i]))
          }
        }
        for (k in 0:(m[l]-1)){
          hpsi3 <- hpsi3+log(1+k*exp(psi[i]))
        }
        v <- c(v,v1+v2)
        hpsi <- c(hpsi,hpsi1+hpsi2-hpsi3)
      }
    }

    svs <- s*v*s
    w <- sqrt(svs)
    #w[which(w<10^(-6))] <- 10^(-6)

    # W <- diag(c(svs))
    # W[which(W<10^(-6))] <- 10^(-6)

    # Compute H (Expected Hessian Matrix)
    # FASTEN CALCULATIONS
    # H1 <- cbind(t(X)%*%W%*%X,t(X)%*%W%*%Z)
    # H2 <- cbind(t(Z)%*%W%*%X,t(Z)%*%W%*%Z+D.)
    H1 <- cbind(crossprod(X*w),crossprod(X*w,Z*w))
    H2 <- cbind(crossprod(Z*w,X*w),crossprod(Z*w)+D.)
    H <- rbind(H1,H2)
    eig <- svd(H)$d
    H.. <- sum(log(eig))
    out <- -(sum(hpsi)-(1/2)*H..)

    return(out)

  }

  # Method for the optimization
  Method <- ifelse(nDim==1,"Brent","L-BFGS-B")
  OLDpsi <- log(OLDphi)
  mle.psi <- optim(OLDpsi,function.psi,hessian=TRUE,method=Method,lower=rep(-10,nDim),upper=rep(4,nDim),control=list(pgtol=1e-3))
  psi <- mle.psi$par

  if (min(psi)==-10){
    psi[order(psi)[1]] <- runif(1,-10,-8)
    #psi[which(psi==-100)] <- -100
    cat("-Phi estimation has not converged reaching the lower bound-\n")
  }
  phi <- exp(psi)

  #-----------------------------------------------------------------------------
  # ESTIMATION OF SIGMA
  #-----------------------------------------------------------------------------
  # Compute v because it is the same for all sigma score equation
  v <- NULL
  for (i in 1:nDim){
    for (l in ((i-1)*nObs+1):(i*nObs)){
      v1 <- 0
      v2 <- 0
      if (y[l]==0){
        v1 <- 0
      }else{
        for (k in 0:(y[l]-1)){
          v1 <- v1+1/(p[l]+k*phi[i])^2
        }
      }
      if (y[l]==m[l]){
        v2 <- 0
      }else{
        for (k in 0:(m[l]-y[l]-1)){
          v2 <- v2+1/(1-p[l]+k*phi[i])^2
        }
      }
      v <- c(v,v1+v2)
    }
  }


  # We get the score equation terms where are the same for all sigma
  # SVS <- diag(s*v*s)
  svs <- sqrt(s*v*s)
  # SVS[which(SVS<10^(-6))] <- 10^(-6) Cambian considerablemente las estimaciones! X
  # FASTEN CALCULATIONS
  # P. <- t(Z)%*%SVS%*%Z-t(Z)%*%SVS%*%X%*%solve(t(X)%*%SVS%*%X)%*%t(X)%*%SVS%*%Z
  MAT1 <- solve(crossprod(X*svs))
  MAT2 <- crossprod(X*svs,Z*svs)
  P. <- crossprod(Z*svs)-crossprod(MAT2,crossprod(MAT1,MAT2))
  #P.[which(P.<10^(-6))] <- 10^(-6) Cambian considerablemente las estimaciones! X


  function.sigma <- function(sigma){

    nRand.previous <- 0
    d. <- NULL
    out <- NULL

    for (i in 1:nComp){
      d. <- c(d.,rep(1/(sigma[i])^2,nRandComp[i]))
    }

    # Compute the trace
    D. <- diag(c(d.))
    P <- solve(P.+D.)
    trace <- diag(P)

    for(i in 1:nComp){

      # SIMPLIFIED
      out <- c(out,-nRandComp[i]*(sigma[i])^2+t(u[seq(nRand.previous+1,nRand.previous+nRandComp[i])])%*%u[seq(nRand.previous+1,nRand.previous+nRandComp[i])]+sum(trace[seq(nRand.previous+1,nRand.previous+nRandComp[i])]))

      # ORIGINAL
      #out <- c(out,-nRandComp[i]/(sigma[i])+t(u[seq(nRand.previous+1,nRand.previous+nRandComp[i])])%*%u[seq(nRand.previous+1,nRand.previous+nRandComp[i])]/(sigma[i])^3+sum(trace[seq(nRand.previous+1,nRand.previous+nRandComp[i])])/(sigma[i])^3)

      nRand.previous <- nRand.previous+nRandComp[i]
    }
    out
  }

  mle.sigma <- multiroot(function.sigma,start=OLDall.sigma,atol=0.01,ctol=0.01,rtol=0.01)
  all.sigma <- mle.sigma$root

  #-----------------------------------------------------------------------------
  # VARIANCES
  #-----------------------------------------------------------------------------
  # Variance of psi
  psi.var <- 1/diag(mle.psi$hessian)
  # Variance of sigma: in the uniroot equation the derivative of the score equation has been simplifyed, therefore we define the function without
  # any simplification to calculate the derivative (2nd derivative of the score equation).
  # We have simplifyed the derivative of the score equation because if not sometimes it gets errors as there are multipli solutions of sigma (0).

  function.sigma <- function(sigma){
    nRand.previous <- 0
    d. <- NULL
    out <- NULL
    for (i in 1:nComp){
      d. <- c(d.,rep(1/(sigma[i])^2,nRandComp[i]))
    }
    # Compute the trace
    D. <- diag(c(d.))
    P <- solve(P.+D.)
    trace <- diag(P)
    for(i in 1:nComp){
      out <- c(out,-nRandComp[i]/(sigma[i])+t(u[seq(nRand.previous+1,nRand.previous+nRandComp[i])])%*%u[seq(nRand.previous+1,nRand.previous+nRandComp[i])]/(sigma[i])^3+sum(trace[seq(nRand.previous+1,nRand.previous+nRandComp[i])])/(sigma[i])^3)
      nRand.previous <- nRand.previous+nRandComp[i]
    }
    out
  }
  sigma.var <- -1/grad(function.sigma,all.sigma)

  #-----------------------------------------------------------------------------
  #OUTPUT
  #-----------------------------------------------------------------------------
  out <- list(phi=phi,all.sigma=all.sigma,psi=psi,psi.var=psi.var,all.sigma.var=sigma.var)
  return(out)
}