Last updated: 2017-01-02

Code version: 55e11cf8f7785ad926b716fb52e4e87b342f38e1

# Pre-requisites

include the most complex concepts required to understand the material.

# Overview

Suppose we have a logistic regression $$Y_i | X_i \sim Bern(p_i)$$ where $log(p_i/(1-p_i)) = \mu + \theta X_i.$

We will assume that $$X_i \in {-1,+1}$$, and assume priors for $$\mu$$ and $$\theta$$: $\mu \sim N(0,100)$ $\theta \sim N(0,1)$

For illustration we simulate data where $$\mu=\theta=0$$:

x = sample(c(-1,1),1000,replace=TRUE)
y = rbinom(1000,1,0.5)

#b is a vector b=(mu,theta)
#loglikelihood for logistic regression
loglik = function(b){
eta = b+b*x
p = exp(eta)/(1+exp(eta))
return(sum(log(y*p+(1-y)*(1-p))))
}

#b is a vector b=(mu,theta)
log_prior = function(b){
return(dnorm(b,0,10, log=TRUE)+dnorm(b,0,1,log=TRUE))
}

#b is a vector b=(mu,theta)
log_post = function(b){
return(loglik(b)+log_prior(b))
}

Let’s compute a 95% CI for $$\theta$$. First try a discrete grid

Note: This is still a work in progress.

m = seq(-10,10,length=100)
t = seq(-2,2,length=100)
df = expand.grid(m=m,t=t)
head(df)
#df = c(df,dplyr::ddply(df,log_post))

# Examples

## Session information

sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.5 LTS

locale:
 LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 LC_PAPER=en_US.UTF-8       LC_NAME=C
 LC_ADDRESS=C               LC_TELEPHONE=C
 LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
 stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 knitr_1.14     MASS_7.3-45    expm_0.999-0   Matrix_1.2-7.1
 rmarkdown_1.1

loaded via a namespace (and not attached):
 Rcpp_0.12.7     lattice_0.20-34 gtools_3.5.0    digest_0.6.9
 assertthat_0.1  mime_0.4        R6_2.1.2        grid_3.3.2
 xtable_1.8-2    formatR_1.4     magrittr_1.5    evaluate_0.9
 stringi_1.1.1   tools_3.3.2     stringr_1.0.0   shiny_0.13.2
 httpuv_1.3.3    yaml_2.1.13     htmltools_0.3.5 tibble_1.2     

This site was created with R Markdown