demo/old/delta_by_hand.R

In harvard-P01/causalMediation: Causal Mediation Analysis

df <- data.frame('smoking'    = c(0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0),
                 'lbw'        = c(0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                 'death'      = c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
                 'drinking'   = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0),
                 'agebelow20' = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
)

data = df
outcome = 'death'
treatment = 'smoking'
mediator = 'lbw'
covariates = c('drinking', 'agebelow20')
interaction = TRUE

mediator.basic <- paste(mediator, treatment, sep=' ~ ')
outcome.basic  <- paste(paste(outcome, treatment, sep=' ~ '), mediator, sep=' + ')

if(interaction==TRUE){
  outcome.basic <- paste(outcome.basic, paste(treatment, mediator, sep = '*'), sep = ' + ')
}

if (length(covariates) == 0){
  mediator.formula <- mediator.basic
  outcome.formula  <- outcome.basic
} else {
  mediator.formula <- paste(mediator.basic, paste(covariates, collapse = " + "), sep = ' + ')
  outcome.formula  <- paste(outcome.basic,  paste(covariates, collapse = " + "), sep = ' + ')
}


### FIXME: hardcode to validate with SAS macro
mediator.binary <- all(unique(data[, mediator]) %in% 0:1)
outcome.binary <- all(unique(data[, outcome])  %in% 0:1)

##----- Delta method
if (! mediator.binary){
  mediator.regression <- lm(mediator.formula, data = data)
}else{
  mediator.regression <- glm(mediator.formula, family = binomial(), data = data)
}

if(!outcome.binary){
  outcome.regression  <- lm(outcome.formula, data = data)
}else{
  outcome.regression  <- glm(outcome.formula, family = binomial(), data = data)
}

## Store coefficients from regression
betas  <- coefficients(mediator.regression)
thetas <- coefficients(outcome.regression)

## Store covariances from regression
vcov_betas <- vcov(mediator.regression)
vcov_thetas <- vcov(outcome.regression)

## Build block diagonal matrix
vcov_block <- bdiag(vcov_betas, vcov_thetas)

cde <- CDE_bin(thetas, treatment, mediator)
nde <- NDE_binbin(betas, thetas, treatment, mediator, covariates)
nie <- NIE_binbin(betas, thetas, treatment, mediator, covariates)
tde <- nde*nie

attach(data)
m <- 1
a <- 1
a_star <- 0


length(thetas)
paste("x")

exp()

g <- exp((thetas[2]+thetas[length(thetas)]*m)*(a-a_star))

deltamethod(~ exp((x2+x6*m)*(a-a_star)), thetas, vcov_thetas)
s <- CDE_bin_delta(thetas, "smoking", "lbw", interaction = FALSE)
s

deltamethod(as.formula(s), thetas, vcov_thetas)

syms <- paste("x", 1:n, sep = "")

deltamethod (~ 3, thetas, vcov_thetas) 
deltamethod (~ 1 / (x1 + x2), estmean, estvar) 

# attr(mediator.regression, "terms")

## Simple linear regression, E(y) = alpha + beta x
set.seed(02138)
x <- 1:100
y <- rnorm(100, 4*x, 5)
toy.lm <- lm(y ~ x)
estmean <- coef(toy.lm)
estvar <- summary(toy.lm)$cov.unscaled * summary(toy.lm)$sigma^2

## Estimate of (1 / (alphahat + betahat))
1 / (estmean[1] + estmean[2])
## Approximate standard error
deltamethod (~ 1 / (x1 + x2), estmean, estvar) 

# Simple linear regression, E(y) = alpha + beta x 
set.seed(02138)
a <- 1:100
b <- rnorm(100, 4*x, 5)
toy.lm <- lm(b ~ a)
estmean <- coef(toy.lm)
estvar <- summary(toy.lm)$cov.unscaled * summary(toy.lm)$sigma^2

## Estimate of (1 / (alphahat + betahat))
1 / (estmean[1] + estmean[2])
## Approximate standard error
deltamethod (~ 1 / (x1 + x2), estmean, estvar)