R/BBmm_EffectsEst_multiroot.R

In idaejin/PROreg: Patient Reported Outcomes Regression Analysis

EffectsEst.multiroot <- function(y,m,beta,u,phi,D.,X,Z){

  # number of observations
  nObs <- length(y)

  # number of fixed effects
  q <- length(beta)

  # number of random effects
  nRand <- length(u)

  # Initial values
  oldbeta <- beta
  oldu <- u

  # Defining the score equation functions
  EffectsEst.function <- function(x){
    beta <- x[1:q]
    u <- x[(q+1):length(x)]
    p <- 1/(1+exp(-(X%*%beta+Z%*%u)))
    for (i in 1:nObs){
      if (p[i]==0){
        p[i] <- 0.001
      }
      if (p[i]==1){
        p[i] <- 0.999
      }
    }
    S <- diag(c(p*(1-p)))

    t <- NULL
    for (j in 1:nObs){
      t1 <- 0
      t2 <- 0
      if (y[j]==0){
        t1 <- 0
      }else{
        for (k in 0:(y[j]-1)){
          t1 <- t1+1/(p[j]+k*phi)
        }
      }
      if (y[j]==m[j]){
        t2 <- 0
      }else{
        for (k in 0:(m[j]-y[j]-1)){
          t2 <- t2+1/(1-p[j]+k*phi)
        }
      }
      t <- c(t,t1-t2)
    }
    c(F1=t(X)%*%S%*%t,F2=t(Z)%*%S%*%t-D.%*%u)
  }

  # Geting the mle
  start <- c(c(oldbeta),c(oldu))
  mle <- multiroot(f=EffectsEst.function,start=start,ctol=0.01,atol=0.01,rtol=0.01,useFortran = TRUE)

  # Calculing the values
  beta <- mle$root[1:q]
  names(beta) <- colnames(X)
  u <- mle$root[(q+1):(nRand+q)]

  # return
  out <- list(fixed.est=beta,random.est=u)

}