R/hess.R

In ui: Uncertainty Intervals and Sensitivity Analysis for Missing Data

#' Hessian for the loglikelihood used by ui_probit 
#'
#' This function derives the hessian in order for \code{\link{ui_probit}} to run faster.
#' @param par Coefficients.
#' @param rho Rho.
#' @param Xz Covariate matrix for missingness.
#' @param Xy Covariate matrix for outcome.
#' @param z Missing or not.
#' @param y Outcome.
#' @import mvtnorm
#' @importFrom stats pnorm dnorm
#' @export

hess<- function(par,rho,Xz = Xz, Xy = Xy, y = y, z = z){
d<-dim(Xy)[2]
beta<-par[1:d]
delta<-par[(d+1):length(par)]

q <- 2*y - 1
w1<-(q*tcrossprod(beta,Xy))[z==1]
dx<-tcrossprod(delta,Xz)

q1<-q[z==1]
dx1<-dx[z==1]
Rhos1<-q1*rho
dx0<-dx[z==0]

n1<-sum(z==1)
Phi2<-vector(length=n1)
for(i in 1:n1){
		Phi2[i]<-pmvnorm(lower=-Inf,upper=c(w1[i],dx1[i]), 
						mean=c(0,0),corr=rbind(c(1,Rhos1[i]),c(Rhos1[i],1)))}
						
pnorm1 <- pnorm((dx1 - Rhos1*w1)/sqrt(1 - rho^2))
dnorm1 <- dnorm((dx1 - Rhos1*w1)/sqrt(1 - rho^2))
pnorm2 <- pnorm((w1 - Rhos1*dx1)/sqrt(1 - rho^2))
dnorm2 <- dnorm((w1 - Rhos1*dx1)/sqrt(1 - rho^2))
dnorm_w <- dnorm(w1)
dnorm_dx <- dnorm(dx1)
lambda_dx0<- dnorm(dx0)/(1-pnorm(dx0))


secder_b<- (-dnorm_w/Phi2)*(w1*pnorm1+Rhos1*dnorm1/sqrt(1-rho^2)+ dnorm_w*pnorm1^2/Phi2)
hess_b <- crossprod(c(secder_b)*Xy[z==1,], Xy[z==1,]) 

secder_cross<-(q1*dnorm_w/Phi2)*(dnorm1/sqrt(1-rho^2)-dnorm_dx*pnorm1*pnorm2/Phi2)
hess_cross<-crossprod(c(secder_cross)*Xz[z==1,], Xy[z==1,])

secder_d0<- dx0*lambda_dx0-lambda_dx0^2
secder_d1<- (-dnorm_dx/Phi2)*(dx1*pnorm2+(Rhos1*dnorm2/sqrt(1-rho^2))+dnorm_dx*pnorm2^2/Phi2)
hess_d<-crossprod(c(secder_d0)*Xz[z==0,], Xz[z==0,])+crossprod(c(secder_d1)*Xz[z==1,], Xz[z==1,])

hess<-rbind(cbind(hess_b,t(hess_cross)),cbind(hess_cross,hess_d))
return(hess)
}