R/Robit.R

Defines functions simRobit fitRobit rZeta rWeights rTau rLL rDelta

Documented in fitRobit rDelta rLL rTau rWeights rZeta simRobit

# 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.