causalreg (version 0.1.0)

cglm: Causal generalized linear model

Description

This function does a search for a causal submodel within the generalized linear model provided.

Usage

cglm(
  formula,
  family,
  data,
  alpha = 0.05,
  pval.approx = FALSE,
  B = 100,
  seed = 1,
  search = c("all", "stepwise"),
  ...
)

Value

A glm object of the selected causal submodel.

Arguments

formula

A formula object.

family

A distributional family object. Currently supported options are: binomial and poisson.

data

A data frame containing the variables in the model.

alpha

Significance level for statistical test

pval.approx

If TRUE, chi-squared approximated p-values are calculated. Default is FALSE, in which case p-values are calculated via bootstrap.

B

Number of bootstrap sample when pval.approx=FALSE.

seed

Seed for generating bootstrap samples.

search

If search="stepwise", a greedy forward stepwise search is conducted. Default is search="all", in which case all possible submodels are considered.

...

Further arguments to be passed to the glm function.

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