R/BBest.R

BBest <- function(y,m,method="MM"){
  
  if (sum(as.integer(m)==m)==length(m)){
  } else{
    stop("m must be integer")
  }
  
  if (min(m)<=0){
    stop("m must be positive")
  }
  
  if (is.numeric(y)){
  } else{
    stop("y must be numeric")
  }
  
  if (sum(y%%1==0)==length(y)){
  } else {
    stop("y must be integer")
  }
  
  if ((length(m)==1) | (length(m)==length(y))){
  } else{
    stop("m must be a number, or a vector of the length of y")
  }
  
  if (min(m-y)<0 | min(y) < 0){
    stop("y must be bounded between 0 and m")
  }
  
  if (is.character(method)){
  } else {
    stop("method must be MM or MLE.")
  }
  
  if ((method=="MM") | (method=="MLE")){
  } else {
    stop("method must be MM or MLE.")
  }
  
  #Number of observations
  nObs <- length(y)
  #Balanced
  if (length(m)==1){
    m. <- rep(m,nObs)
    balanced <- "yes"
  } else{
    m. <- m
    if (sum(m[1]==m)==length(m)){
      balanced <- "yes"
    } else{
      balanced <- "no"
    }
  }
  
  #Method of moments
  if (method=="MM"){
    
    p <- sum(y)/sum(m.)
    psi.est <- function(psi){
      out <- (1/(nObs-1))*sum(((y-m.*p)^2)/(m.*p*(1-p)*(1+(m.-1)*(exp(psi)/(1+exp(psi))))))-1
      out
    }
    psi <- uniroot(psi.est,lower = -5,upper =3)$root
    phi <- exp(psi)
  } else{
    
    p <- 0.5
    phi <- 0.5
    psi <- log(phi)
    param <- c(p,phi)
    old.param <- c(Inf,Inf)
    iter <- 0
    
    while(max(abs(param-old.param))>0.001){
      
      old.param <- param
      
      partial.p <- function(p){
        out <- 0
        
        for (j in 1:nObs){
          out1 <- 0
          out2 <- 0
          
          if (y[j]==0){
            out1 <- 0
          }else{
            for (k in 0:(y[j]-1)){
              out1 <- out1+1/(p+k*phi)
            }
          }
          if (y[j]==m.[j]){
            out2 <- 0
          }else{
            for (k in 0:(m.[j]-y[j]-1)){
              out2 <- out2+1/(1-p+k*phi)
            }
          }
          out <- out + (out1-out2)
        }
        return(out)
      }
      
      p <- uniroot(partial.p,lower = 0.001,upper = 0.999)$root
      
      
      partial.psi <- function(psi){
        out <- 0
        for (j in 1:nObs){
          out1 <- 0
          out2 <- 0
          out3 <- 0
          
          if (y[j]==0){
            out1 <- 0
          }else{
            for (k in 0:(y[j]-1)){
              out1 <- out1+(k*exp(psi))/(p+k*exp(psi))
            }
          }
          if (y[j]==m.[j]){
            out2 <- 0
          }else{
            for (k in 0:(m.[j]-y[j]-1)){
              out2 <- out2+(k*exp(psi))/(1-p+k*exp(psi))
            }
          }
          for (k in 0:(m.[j]-1)){
            out3 <- out3 + (k*exp(psi))/(1+k*exp(psi))
          }
          out <- out + (out1+out2-out3)
        }
        return(out)
      }
      
      psi <- uniroot(partial.psi,lower = -5,upper = 3)$root
      
      phi <- exp(psi)
      param <- c(p,phi)
      
      iter <- iter+1
    }
    p <- param[1]
    phi <- param[2]
  }
  
  
  # The variance is estimated through the Fisher Information Matrix:
  partial.p <- function(p){
    out <- 0
    
    for (j in 1:nObs){
      out1 <- 0
      out2 <- 0
      
      if (y[j]==0){
        out1 <- 0
      }else{
        for (k in 0:(y[j]-1)){
          out1 <- out1+1/(p+k*phi)
        }
      }
      if (y[j]==m.[j]){
        out2 <- 0
      }else{
        for (k in 0:(m.[j]-y[j]-1)){
          out2 <- out2+1/(1-p+k*phi)
        }
      }
      out <- out + (out1-out2)
    }
    return(out)
  }
  
  partial.psi <- function(psi){
    out <- 0
    for (j in 1:nObs){
      out1 <- 0
      out2 <- 0
      out3 <- 0
      
      if (y[j]==0){
        out1 <- 0
      }else{
        for (k in 0:(y[j]-1)){
          out1 <- out1+(k*exp(psi))/(p+k*exp(psi))
        }
      }
      if (y[j]==m.[j]){
        out2 <- 0
      }else{
        for (k in 0:(m.[j]-y[j]-1)){
          out2 <- out2+(k*exp(psi))/(1-p+k*exp(psi))
        }
      }
      for (k in 0:(m.[j]-1)){
        out3 <- out3 + (k*exp(psi))/(1+k*exp(psi))
      }
      out <- out + (out1+out2-out3)
    }
    return(out)
  }
  
  var.p <- -1/grad(partial.p,p)
  var.psi <- -1/grad(partial.psi,psi)
  
  out <- list(p=p, phi=phi,
              pVar=var.p,psi=psi,psiVar=var.psi,
              m=m,balanced=balanced,method=method)
  
  class(out) <- "BBest"
  
  out$call <- match.call()
  
  out
}


print.BBest <- function(x,...){
  
  cat("The probability parameter of the beta-binomial distribution:",x$p,"\n")
  cat("The dispersion parameter of the beta-binomial distribution:",x$phi,"\n")
  if (length(x$m)==1){
    cat("\nBalanced data, maximum score number:",x$m,"\n")
  } else {
    cat("\nNo balanced data.\n")
  }
}


summary.BBest <- function(object,...){
  p.se <- sqrt(object$pVar)
  p.coef <- cbind(object$p,p.se)
  colnames(p.coef) <- c("Estimate","StdErr")
  
  psi.se <- sqrt(object$psiVar)
  psi.coef <- cbind(object$psi,psi.se)
  colnames(psi.coef) <- c("Estimate","StdErr")
  
  
  res <- list(call=object$call,coefficients=cbind(p=object$p,phi=object$phi),
              p.coefficients=p.coef, psi.coefficients=psi.coef,
              m=object$m)
  
  class(res) <- "summary.BBest"
  res
}


print.summary.BBest <- function(x,...){
  cat("Call:\t")
  print(x$call)
  cat("\n")
  cat("Probability parameter estimation:\n")
  print(x$p.coefficients)
  cat("\n")
  cat("Log dispersion parameter estimation:\n")
  print(x$psi.coefficients)
  
  if (length(x$m)==1){
    cat("\nBalanced data, maximum score number:",x$m,"\n")
  } else {
    cat("\nNo balanced data.\n")
  }
}

Try the HRQoL package in your browser

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

HRQoL documentation built on May 2, 2019, 5:42 a.m.