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

#### Documented in fitRobitrDeltarLLrTaurWeightsrZetasimRobit

# Purpose : Estimation of Robit regression model

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

#' Delta Matrix
#'
#' Delta matrix for robit IRLS
#' @param eta Linear predictor
#' @param nu Degrees of freedom
#'
#' @importFrom stats dt

rDelta = function(eta,nu=7){
a = dt(x=eta,df=nu);
return(1/a);
}

#' Likelihood
#'
#' Robit model log likelihood
#' @param y Outcome
#' @param eta Linear predictor
#' @param nu Degrees of freedom
#'
#' @importFrom stats pt

rLL = function(y,eta,nu=7){
# Ones
l1 = sum(pt(q=eta[y==1],df=nu,log.p=T));
# Zeros
l0 = sum(pt(q=-1*eta[y==0],df=nu,log.p=T));
return(l1+l0);
}

#' Tau1
#'
#' Expectation of \eqn{\tau} given y in the robit model
#' @param y Outcome
#' @param eta Linear predictor
#' @param nu Degrees of freedom
#' @importFrom stats pt

rTau = function(y,eta,nu=7){
n = length(y);
# Numerator
a1 = -1*sqrt(1+2/nu)*eta;
a = y-(2*y-1)*pt(q=a1,df=(nu+2));
# Denominator
b = y-(2*y-1)*pt(q=-1*eta,df=nu);
# Format
Out = diag(as.numeric(a/b),ncol=n);
return(Out);
}

#' Weight Matrix
#'
#' Weight matrix for robit IRLS
#' @param eta Linear predictor
#' @param nu Degrees of freedom
#'
#' @importFrom stats dt pt

rWeights = function(eta,nu=7){
n = length(eta);
a = as.numeric(dt(x=eta,df=nu)^2);
b = as.numeric(pt(q=eta,df=nu)*pt(q=-1*eta,df=nu));
return(diag(a/b,nrow=n));
}

#' Robit Zeta
#'
#' Expectation of z given y in the robit model
#' @param y Outcome
#' @param eta Linear predictor
#' @param nu Degrees of freedom
#' @importFrom stats dt pt

rZeta = function(y,eta,nu){
# Numerator
a = (2*y-1)*dt(x=eta,df=nu);
# Denominator
b1 = -1*sqrt(1+2/nu)*eta;
b = y-(2*y-1)*pt(q=b1,df=(nu+2));
return(eta+(a/b));
}

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

#' Robit 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
#' @param nu Degrees of freedom
#'
#' @importFrom stats model.matrix pnorm
#' @export

fitRobit = function(y,X,b0,m1=100,m2=10,nu=7){
# 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 = rLL(y=y,eta=eta0);
for(i in 1:m1){
# Expectation of tau given y
T1 = rTau(y=y,eta=eta0,nu=nu);
# Expectation of z given y
z1 = rZeta(y=y,eta=eta0,nu=nu);
# Update expanded beta
b1x = wls(z=z1,X=Dx,W=T1);
# Update sigma
s1 = rSigma(nu=nu,z1=z1,eta=eta0,T1=T1,X=Dx);
# Update alpha
a1 = sum(diag(T1))/n;
# Update reduced beta
b1 = sqrt(a1)*b1x/s1;
# Check for likelihood improvement
eta1 = Dx %*% b1;
ll1 = rLL(y=y,eta=eta1,nu=nu);
if(ll1<=ll0){break};
# Update
ll0 = ll1;
b0 = b1;
eta0 = eta1;
}
# IRLS
for(j in 1:m2){
# Delta
Delta = rDelta(eta=eta0,nu=nu);
# Pi
Pi = pt(eta0,df=nu);
# Working vector
z = eta0 + Delta*(y-Pi);
# Weights
W = rWeights(eta=eta0,nu=nu);
# WLS update
b1 = wls(z=z,X=Dx,W=W);
# New linear predictor
eta1 = Dx %*% b1;
# Check for likelihood improvement
ll1 = rLL(y=y,eta=eta1,nu=nu);
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
########################

#' Robit Simulation
#'
#' Simulate binary outcomes from a probit model with linear predictor \eqn{\eta}.
#' @param eta Linear predictor
#' @param nu Degrees of freedom
#'
#' @importFrom stats rgamma rnorm
#' @export

simRobit = function(eta,nu=7){
n = length(eta);
tau = rgamma(n=n,shape=nu/2,rate=nu/2);
z = rnorm(n=n,mean=eta,sd=1/tau);
y = 1*(z>0);
return(y);
}

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