R/cgam.R

In causalreg: Causal Generalized Linear Models

#' Causal generalized additive model
#'
#' This function does a search for a causal submodel within the generalized additive 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 gam function.
#' @return A gam object of the selected causal submodel.
#' @importFrom stats BIC as.formula coef glm pchisq reformulate residuals terms.formula update.formula
#' @importFrom utils combn
#' @importFrom mgcv gam
#' @examples
#' ##############################
#' #causal Poisson gam##########
#' n<-1000
#' set.seed(123)
#' X1<-rnorm(n,0,1)
#' Y<-rpois(n,exp(sin(X1)))
#' X2<-log(Y+1)+rnorm(n,0,0.5)
#' data<-data.frame(cbind(X1, X2, Y))
#' cm_all<-cgam(Y ~ s(X1)+s(X2),"poisson",data,pval.approx=TRUE,search="all")
#' cm_all$model.opt
#' cm_step<-cgam(Y ~ s(X1)+s(X2),"poisson",data,pval.approx=TRUE,search="stepwise")
#' cm_step$model.opt
#' \donttest{
#' #bigger simulation with 7 covariates
#' set.seed(123)
#' n<-1000
#' X1<-rnorm(n=n,sd=sqrt(0.04))
#' X2<-X1+rnorm(n=n,sd=sqrt(0.04))
#' X3<-X1+X2+rnorm(n=n,sd=sqrt(0.04))
#' m<-sin(X2*5)+X3^3
#' Z<-m+rnorm(n=n,sd=sqrt(0.04))
#' X4<-X2+rnorm(n=n,sd=sqrt(0.04))
#' X5<-Z+rnorm(n=n,sd=sqrt(0.04))
#' X6<-Z+rnorm(n=n,sd=sqrt(0.04))
#' X7<-X6+rnorm(n=n,sd=sqrt(0.04))
#' Y<-qpois(pnorm(Z, mean = m, sd = sqrt(0.04)), lambda=exp(m))
#' dat<-data.frame(cbind(X1, X2, X3, X4, X5, X6, X7,Y))
#' fml<- Y~s(X1)+s(X2)+s(X3)+s(X4)+s(X5)+s(X6)+s(X7)
#' mod.all <-cgam(fml,"poisson",dat,pval.approx=TRUE,search="all")
#' mod.all$model.opt
#' mod.step <-cgam(fml,"poisson",dat,pval.approx=TRUE,search="stepwise")
#' mod.step$model.opt
#' ####################################
#' #causal logistic gam################
#' n<-1000
#' 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<-cgam(Y ~ s(X1)+s(X2),"binomial",data,pval.approx=FALSE,search="all")
#' cm_all$model.opt
#' cm_step<-cgam(Y ~ s(X1)+s(X2),"binomial",data,pval.approx=FALSE,search="stepwise")
#' cm_step$model.opt
#' }
#' @export
cgam<-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<-cgam_all(formula, family, data, alpha,pval.approx,B,...)
  }
  else
  {
    res<-cgam_step(formula, family, data, alpha,pval.approx,B,...)
  }
  return(res)
}