R/grr.R

In ui: Uncertainty Intervals and Sensitivity Analysis for Missing Data

#' Gradient for the loglikelihood used by ui_probit 
#'
#' This function derives the gradient 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
grr<- 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
bx<-tcrossprod(beta,Xy)
bx1<-bx[z==1]
w1<-q[z==1]*bx1
dx<-tcrossprod(delta,Xz)
q1<-q[z==1]
Rhos1<-q1*rho
dx1<-dx[z==1]

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)))}

gr_b<-q1*dnorm(w1)*pnorm((dx1-rho*bx1)/sqrt(1-rho^2))/Phi2

gr_beta<-crossprod(gr_b,Xy[z==1,])

gr_d1 <- dnorm(dx1)*pnorm(q1*(bx1-rho*dx1)/sqrt(1-rho^2))/Phi2 
gr_d0 <- -dnorm(dx[z==0])/(1-pnorm(dx[z==0]))

gr_delta<-crossprod(gr_d1,Xz[z==1,])+crossprod(gr_d0,Xz[z==0,])

return(c(gr_beta,gr_delta))
}