R/cglm.R

In causalreg: Causal Generalized Linear Models

#' Causal generalized linear model
#'
#' This function does a search for a causal submodel within the generalized linear model provided.
#'
#' @param formula A formula object.
#' @param family A distributional family object. Currently supported options are: binomial and poisson.
#' @param data A data frame containing the variables in the model.
#' @param alpha Significance level for statistical test
#' @param pval.approx If TRUE, chi-squared approximated p-values are calculated. Default is FALSE, in which case p-values are calculated via bootstrap.
#' @param B Number of bootstrap sample when pval.approx=FALSE.
#' @param seed Seed for generating bootstrap samples.
#' @param search If search="stepwise", a greedy forward stepwise search is conducted. Default is search="all", in which case all possible submodels are considered.
#' @param ... Further arguments to be passed to the glm function.
#' @return A glm object of the selected causal submodel.
#' @importFrom stats BIC as.formula coef glm pchisq reformulate residuals terms.formula update.formula
#' @importFrom utils combn
#' @examples
#' ###################################
#' #causal Poisson glm#################
#' n<-1000
#' set.seed(123)
#' X1<-rnorm(n,0,1)
#' Y<-rpois(n,exp(X1))
#' X2<-log(Y+1)+rnorm(n,0,0.3)
#' data<-data.frame(cbind(X1, X2, Y))
#' cm_all<-cglm(Y ~ X1+X2,"poisson",data,pval.approx=TRUE,search="all")
#' cm_all$model.opt
#' cm_step<-cglm(Y ~ X1+X2,"poisson",data,pval.approx=TRUE,search="stepwise")
#' cm_step$model.opt
#' \donttest{
#' ##########################
#' #causal logistic glm#######
#' n<-2000
#' set.seed(123)
#' X1<-rnorm(n,0,1)
#' Y<-rbinom(n,1,exp(X1)/(1+exp(X1)))
#' flip<-rbinom(n,1,0.1)
#' X2<-(1-flip)*Y+rnorm(n,0,0.3)
#' data<-data.frame(cbind(X1, X2, Y))
#' cm_all<-cglm(Y ~ X1+X2,"binomial",data,pval.approx=FALSE,search="all")
#' cm_all$model.opt
#' cm_step<-cglm(Y ~ X1+X2,"binomial",data,pval.approx=FALSE,search="stepwise")
#' cm_step$model.opt
#' #bigger simulation with 5 covariates
#' set.seed(12)
#' n<-3000
#' X1<-rnorm(n,0,1)
#' X2<-rnorm(n,X1,0.5)
#' X3<-rnorm(n,0,1)
#' X4<-rnorm(n,X2,.5)
#' Y<-rbinom(n,1,exp(.8*X2-.9*X3)/(1+exp(.8*X2-.9*X3)))
#' flip<-rbinom(n,1,0.1)
#' X5<-(1-flip)*Y+flip*(1-Y)+rnorm(n,0,.3)
#' dat<-data.frame(cbind(X1, X2, X3, X4, X5,Y))
#' mod.all <-cglm(Y~X1+X2+X3+X4+X5,"binomial",dat,pval.approx=FALSE,search="all")
#' mod.all$model.opt
#' mod.step <-cglm(Y~X1+X2+X3+X4+X5,"binomial",dat,pval.approx=FALSE,search="stepwise")
#' mod.step$model.opt
#' }
#' @export
cglm<-function(formula, family, data, alpha=0.05,pval.approx=FALSE,B=100,seed=1,search=c("all","stepwise"),...)
{
  set.seed(seed)

  if(search[1]=="all")
  {
    res<-cglm_all(formula, family, data, alpha,pval.approx,B,...)
  }
  else
  {
    res<-cglm_step(formula, family, data, alpha,pval.approx,B,...)
  }
  return(res)
}