# R/Probit.R In zrmacc/BinaryModels: Binary Regression Models

#### Documented in fitProbitpDeltapLLpWeightspZeta1pZeta2simProbittwoNorm

#' @useDynLib BinaryModels, .registration = TRUE
#' @importFrom Rcpp sourceCpp
NULL

########################
# Auxiliary Functions
########################

#' Delta Matrix
#'
#' Delta matrix for probit IRLS
#' @param eta Linear predictor
#'
#' @importFrom stats dnorm

pDelta = function(eta){
a = dnorm(eta);
return(1/a);
}

#' L2 Norm
#'
#' @param b Vector

twoNorm = function(b){
a = sum(b^2);
return(sqrt(a));
}

#' Likelihood
#'
#' Probit model log likelihood
#' @param y Outcome
#' @param eta Linear predictor
#'
#' @importFrom stats pnorm

pLL = function(y,eta){
# Ones
l1 = sum(pnorm(eta[y==1],log.p=T));
# Zeros
l0 = sum(pnorm(-eta[y==0],log.p=T));
return(l1+l0);
}

#' Weight Matrix
#'
#' Weight matrix for probit IRLS
#' @param eta Linear predictor
#'
#' @importFrom stats dnorm pnorm

pWeights = function(eta){
n = length(eta);
a = as.numeric(dnorm(eta)^2);
b = as.numeric(pnorm(eta)*pnorm(-eta));
return(diag(a/b,nrow=n));
}

#' Probit Zeta1
#'
#' Expectation of z given y in the probit model
#' @param y Outcome
#' @param eta Linear predictor
#' @param alpha Expansion parameter
#'
#' @importFrom stats dnorm pnorm

pZeta1 = function(y,eta,alpha){
n = length(y);
k = rep(0,n);
# Ones
k[y==1] = alpha*dnorm(eta[y==1]/alpha)/pnorm(eta[y==1]/alpha);
# Zeros
k[y==0] = -1*alpha*dnorm(eta[y==0]/alpha)/pnorm(-1*eta[y==0]/alpha);
return(eta+k);
}

#' Zeta2
#'
#' Expectation of \eqn{z^2} given y in the probit model
#' @param y Outcome
#' @param eta Linear predictor
#' @param alpha Expansion parameter
#'
#' @importFrom stats dnorm pnorm

pZeta2 = function(y,eta,alpha){
n = length(y);
k = rep(0,n);
k[y==1] = eta[y==1]*alpha*dnorm(eta[y==1]/alpha)/pnorm(eta[y==1]/alpha);
k[y==0] = -1*eta[y==0]*alpha*dnorm(eta[y==0]/alpha)/pnorm(-1*eta[y==0]/alpha);
a = (eta^2+alpha^2)
return(a+k);
}

########################
# Model Fitting
########################

#' Probit Regression
#'
#' @param y Binary outcome
#' @param X Design matrix
#' @param b0 Initial estimate of \eqn{\beta}, defaults to zero.
#' @param m1 Maximum iterations of PX-EM
#' @param m2 Maximum iterations of IRLS
#'
#' @importFrom stats model.matrix pnorm
#' @export

fitProbit = function(y,X,b0,m1=100,m2=10){
# Convert X to model matrix
Dx = model.matrix(~.,data=data.frame(X));
k = ncol(Dx);
# Check for initial estimate of beta
if(missing(b0)){b0 = rep(0,k)};
# Initialization
eta0 = Dx %*% b0;
ll00 = ll0 = pLL(y=y,eta=eta0);
for(i in 1:m1){
# Expectation of z given y
z1 = pZeta1(y=y,eta=eta0,alpha=1);
# Update expanded beta
b1x = pBeta(z1=z1,X=Dx);
# Expectation of z2 given y
z2 = pZeta2(y=y,eta=eta0,alpha=1);
# Update alpha
a1 = pAlpha(z1=z1,z2=z2,eta=eta0);
# Update reduced beta
b1 = b1x/a1;
# Check for likelihood improvement
eta1 = Dx %*% b1;
ll1 = pLL(y=y,eta=eta1);
if(ll1<=ll0){break};
# Update
ll0 = ll1;
b0 = b1;
eta0 = eta1;
}
# IRLS
for(j in 1:m2){
# Delta
Delta = pDelta(eta=eta0);
# Pi
Pi = pnorm(eta0);
# Working vector
z = eta0 + Delta*(y-Pi);
# Weights
W = pWeights(eta=eta0);
# WLS update
b1 = wls(z=z,X=Dx,W=W);
# New linear predictor
eta1 = Dx %*% b1;
# Check for likelihood improvement
ll1 = pLL(y=y,eta=eta1);
if(ll1<=ll0){break};
# Update
ll0 = ll1;
eta0 = eta1;
b0 = b1;
}
# Results
A = Info(X=Dx,W=W);
V = fastInv(A=A);
Out = list("Beta"=b0,"Var"=V);
return(Out);
}

########################
# Simulation Function
########################

#' Probit Simulation
#'
#' Simulate binary outcomes from a probit model with linear predictor \eqn{\eta}.
#' @param eta Linear predictor
#'
#' @importFrom stats rnorm
#' @export

simProbit = function(eta){
n = length(eta);
z = rnorm(n=n,mean=eta,sd=1);
y = 1*(z>0);
return(y);
}

zrmacc/BinaryModels documentation built on Jan. 14, 2018, 9:09 a.m.