R/ui_probit.R

In ui: Uncertainty Intervals and Sensitivity Analysis for Missing Data

#' Uncertainty intervals for probit regression
#'
#' This function allows you to derive uncertainty intervals for probit regression
#'	when there is missing data in the binary outcome. The uncertainty intervals
#'	can be used as a sensitivity analysis to ignorability (missing at random), and
#'	are derived by maximum likelihood. Note that \code{rho}=0 render the same results as
#'	a complete case analysis.
#' @param out_formula Formula for outcome regression.
#' @param mis_formula Formula for missingness mechanism. If NULL the same covariates as in the outcome regression will be used.
#' @param data data.frame containing the variables in the formula.
#' @param rho Vector containing the values of \code{rho} for which we want to fit the
#'	likelihood.
#' @param progress If TRUE prints out process time for each maximization of the likelihood.
#' @param max_grid Maximum distance between two elements in \code{rho}, if two wide there can
#'	difficulties with convergence of the maximum likelihood.
#' @param alpha Default 0.05 corresponding to a confidence level of 95 for CI and UI.
#' @param method Maximization method to be passed through \code{maxLik}
#'
#'@details In order to visualize the results, you can use \code{\link{plot.uiprobit}}
#'or \code{\link{profile.uiprobit}}.
#'
#' @importFrom  Matrix rankMatrix
#' @return A list containing:
#' \item{coef}{Estimated coefficients (outcome regression) for different values of \code{rho}.}
#' \item{rho}{The values of \code{rho} for which the likelihood is maximized.}
#' \item{vcov}{Covariance matrix.}
#' \item{ci}{Confidence intervals for different values of \code{rho}.}
#' \item{ui}{Uncertainty intervals.}
#' \item{out_model}{Outcome regression model when rho=0.}
#' \item{mis_model}{Regression model for missingness mechanism (selection).}
#' \item{se}{Standard errors from outcome regression.}
#' \item{value}{Value of maximum likelihood for different values of \code{rho}.}
#' \item{y}{Outcome vector.}
#' \item{z}{Indicator variable of observed outcome.}
#' \item{Xy}{Covariate matrix for outcome regression.}
#' \item{Xz}{Covariate matrix for missingness mechanism (selection regression model).}
#' \item{max.info}{Information about the maximization procedure. Includes whether it \code{converged}, \code{message}, \code{method} and \code{number of iterations}.}
#'@author Minna Genbäck
#'@references Genbäck, M., Ng, N., Stanghellini, E., de Luna, X. (2018). Predictors of Decline in Self-reported Health: Addressing Non-ignorable Dropout in Longitudinal Studies of Aging. \emph{European journal of ageing}, 15(2), 211-220.
#' @examples
#'library(MASS)
#'n<-500
#'
#'delta<-c(0.5,0.6,0.1,-1,1)
#'beta<-c(-0.3,-0.5,0,-0.4,-0.3)
#'
#'X<-cbind(rep(1,n),rnorm(n),runif(n),rbinom(n,2,0.5),rbinom(n,1,0.5))
#'x<-X[,-1]
#'rho=0.4
#'error<-mvrnorm(n,c(0,0),matrix(c(1,rho,rho,1),2))
#'
#'zstar<-X%*%delta+error[,1]
#'z<-as.numeric(zstar>0)
#'
#'ystar<-X%*%beta+error[,2]
#'y<-as.integer(ystar>0)
#'y[z==0]<-NA
#'data=data.frame(y=y,x1=x[,1],x2=x[,2],x3=x[,3],x4=x[,4])
#'
#'m<-ui_probit(y~x1+x2+x3+x4,data=data,rho=c(0,0.4))
#'m
#'plot(m)
#'profile(m)
#' @importFrom stats qnorm
#' @export

ui_probit<-function(out_formula,mis_formula=NULL,data,rho=c(-0.3,0.3),progress=TRUE,max_grid=0.1,alpha=0.05,method='NR'){

  cll <- match.call()
   #Warnings
if(!inherits(data,"data.frame")){
	stop('Data must be a data frame')
    }
   rho<-sort(rho)
   if(rho[1]< (-1)|rho[2]>1){
	stop("The limits for rho should be between -1 and 1.")
    }
  # Check to make sure 0 is part of Rho
  if(min(rho)>0 | max(rho)<0){
  rho<-c(0,rho)
  warning('Note that, the rho interval did not cover 0 as recommended. Hence 0 is now added into rho.')
  }
  if(sum(rho==0)==0){
  	rho<-c(rho,0)
  }
  nrho<-length(rho)


  #Make sure that the grid of rho is as fine as specified by max_grid
  for(i in 1:(nrho-1)){
  nri<-ceiling((rho[i+1]-rho[i])/max_grid)
  if(nri>1){
  rho<-c(rho,seq(from=rho[i],to=rho[i+1],by=((rho[i+1]-rho[i])/nri))[2:nri])
  }}

  rho<-sort(rho)
  nrho<-length(rho)

  ml_object<-ml_probit(out_formula,mis_formula,data,rho,progress,method)

  p<-dim(ml_object$coef)[1]-1
  #derive standard errors from vcov matrix
  se<-lapply(ml_object$vcov,function(x) sqrt(diag(x)))
  ml_object$se<-matrix(unlist(se),ncol=nrho)[1:(p+1),]

  zalpha<-qnorm(1-alpha/2)
  ml_object$zalpha<-zalpha

  ml_object$ci<-array(dim=c((p+1),nrho,2),dimnames=list(rownames(ml_object$coef),rho,c('lower','upper')))
  ml_object$ci[,,1]<-ml_object$coef-zalpha*ml_object$se
  ml_object$ci[,,2]<-ml_object$coef+zalpha*ml_object$se
  ml_object$ui<-cbind(apply(ml_object$ci[,,1],1,min),apply(ml_object$ci[,,2],1,max))
  colnames(ml_object$ui)<-c('lower','upper')

  ml_object$call <- cll

  class(ml_object)<-"uiprobit"
  return(ml_object)
}