R/Probit.R

Defines functions simProbit fitProbit pZeta2 pZeta1 pWeights pLL twoNorm pDelta

Documented in fitProbit pDelta pLL pWeights pZeta1 pZeta2 simProbit twoNorm

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