Christian Bartels
Allschwil, June-2014
Efficient generic integration algorithm to determine confidence intervals and p-values for hypothesis testing. Christian Bartels. figshare.
http://dx.doi.org/10.6084/m9.figshare.1054694
# clean up
rm(list = ls())
# print out version information, path, date
R.Version()
## $platform
## [1] "x86_64-w64-mingw32"
##
## $arch
## [1] "x86_64"
##
## $os
## [1] "mingw32"
##
## $system
## [1] "x86_64, mingw32"
##
## $status
## [1] ""
##
## $major
## [1] "3"
##
## $minor
## [1] "1.0"
##
## $year
## [1] "2014"
##
## $month
## [1] "04"
##
## $day
## [1] "10"
##
## $`svn rev`
## [1] "65387"
##
## $language
## [1] "R"
##
## $version.string
## [1] "R version 3.1.0 (2014-04-10)"
##
## $nickname
## [1] "Spring Dance"
date()
## [1] "Thu Jun 12 19:39:53 2014"
Sys.info()
## sysname release version nodename
## "Windows" "7 x64" "build 9200" "IMPERIALDRAMON"
## machine login user effective_user
## "x86-64" "Christian" "Christian" "Christian"
# get libraries
library(rstan)
## Loading required package: Rcpp
## Loading required package: inline
##
## Attaching package: 'inline'
##
## Das folgende Objekt ist maskiert from 'package:Rcpp':
##
## registerPlugin
##
## rstan (Version 2.2.0, packaged: 2014-02-14 04:29:17 UTC, GitRev: 52d7b230aaa0)
library(ggplot2)
library(data.table)
set_cppo("fast")
## mode fast without NDEBUG for compiling C++ code has been set already
# my likelihood function x ... observations or samples of observations col:
# repeated measures (evaluated with the same parameter) row: different
# samples (evaluate with different parameters) par ... data.frame of
# parameters col: different parameters, e.g., mean and sd row: samples
# return vector of log-likelihood of samples
my.ll <- function(x, par) {
return(rowSums(dnorm(x, par$mean, par$sd, log = TRUE)))
}
# my maximum likelihood estimator
my.maxll <- function(x) {
df1 <- do.call(rbind, apply(x, 1, function(x1) {
n <- length(x1)
dft <- data.frame(mean = mean(x1), sd = sd(x1) * sqrt((n - 1)/n)) # max. likelihood parameters (biased)
dft$mll <- my.ll(matrix(x1, nrow = 1), dft) # max. likelihood
return(dft)
}))
return(df1)
}
# my log prior par ... data.frame of parameters col: different parameters,
# e.g., mean and sd row: samples return vector of log-priors of samples
my.prior <- function(par) {
return(-log(par$sd))
}
# my posterior
my.post <- function(x, par) return(my.ll(x, par) + my.prior(par))
# my simulation function par ... data.frame of parameters col: different
# parameters, e.g., mean and sd row: samples n ... number of repeated
# measures
my.sim <- function(n, par) {
m <- nrow(par)
return(matrix(rnorm(n * m, par$mean, par$sd), ncol = n))
}
Data is simulated from a normal distribution.
# make reproducible
set.seed(1)
# generate observed data
n = 10 # number of repeated observations
p0 <- data.frame(mean = 1, sd = 1) # 'unknown truth', not used except for simulation on next line
x <- my.sim(n, p0) # observed data
# list model
modelSample = "mcmc4.model" # model to sample parameters
cbind(system(paste("cat ", modelSample, sep = ""), intern = TRUE))
## [,1]
## [1,] "data {"
## [2,] " int<lower=0> N;"
## [3,] " vector[N] x;"
## [4,] "} "
## [5,] "parameters {"
## [6,] " real m;"
## [7,] " real<lower=0> s;"
## [8,] "} "
## [9,] "model {"
## [10,] " x ~ normal(m, s); // fit"
## [11,] " increment_log_prob(-log(s)); // prior"
## [12,] "} "
# compile model
parSampler <- stan_model(modelSample, verbose = FALSE)
##
## TRANSLATING MODEL 'mcmc4' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'mcmc4' NOW.
## cygwin warning:
## MS-DOS style path detected: C:/PROGRA~1/R/R-31~1.0/etc/x64/Makeconf
## Preferred POSIX equivalent is: /cygdrive/c/PROGRA~1/R/R-31~1.0/etc/x64/Makeconf
## CYGWIN environment variable option "nodosfilewarning" turns off this warning.
## Consult the user's guide for more details about POSIX paths:
## http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
## C:/Users/Christian/Documents/R/win-library/3.1/rstan/include//stansrc/stan/agrad/rev/var_stack.hpp:49:17: warning: 'void stan::agrad::free_memory()' defined but not used [-Wunused-function]
## C:/Users/Christian/Documents/R/win-library/3.1/rstan/include//stansrc/stan/agrad/rev/chainable.hpp:87:17: warning: 'void stan::agrad::set_zero_all_adjoints()' defined but not used [-Wunused-function]
# call Stan
dat <- list(N = n, x = c(x))
parSample <- sampling(parSampler, data = dat, init = 0, iter = 8000, thin = 4,
seed = 123)
## SAMPLING FOR MODEL 'mcmc4' NOW (CHAIN 1).
##
Iteration: 1 / 8000 [ 0%] (Warmup)
Iteration: 800 / 8000 [ 10%] (Warmup)
Iteration: 1600 / 8000 [ 20%] (Warmup)
Iteration: 2400 / 8000 [ 30%] (Warmup)
Iteration: 3200 / 8000 [ 40%] (Warmup)
Iteration: 4000 / 8000 [ 50%] (Warmup)
Iteration: 4800 / 8000 [ 60%] (Sampling)
Iteration: 5600 / 8000 [ 70%] (Sampling)
Iteration: 6400 / 8000 [ 80%] (Sampling)
Iteration: 7200 / 8000 [ 90%] (Sampling)
Iteration: 8000 / 8000 [100%] (Sampling)
## Elapsed Time: 0.044 seconds (Warm-up)
## 0.055 seconds (Sampling)
## 0.099 seconds (Total)
##
## SAMPLING FOR MODEL 'mcmc4' NOW (CHAIN 2).
##
Iteration: 1 / 8000 [ 0%] (Warmup)
Iteration: 800 / 8000 [ 10%] (Warmup)
Iteration: 1600 / 8000 [ 20%] (Warmup)
Iteration: 2400 / 8000 [ 30%] (Warmup)
Iteration: 3200 / 8000 [ 40%] (Warmup)
Iteration: 4000 / 8000 [ 50%] (Warmup)
Iteration: 4800 / 8000 [ 60%] (Sampling)
Iteration: 5600 / 8000 [ 70%] (Sampling)
Iteration: 6400 / 8000 [ 80%] (Sampling)
Iteration: 7200 / 8000 [ 90%] (Sampling)
Iteration: 8000 / 8000 [100%] (Sampling)
## Elapsed Time: 0.044 seconds (Warm-up)
## 0.044 seconds (Sampling)
## 0.088 seconds (Total)
##
## SAMPLING FOR MODEL 'mcmc4' NOW (CHAIN 3).
##
Iteration: 1 / 8000 [ 0%] (Warmup)
Iteration: 800 / 8000 [ 10%] (Warmup)
Iteration: 1600 / 8000 [ 20%] (Warmup)
Iteration: 2400 / 8000 [ 30%] (Warmup)
Iteration: 3200 / 8000 [ 40%] (Warmup)
Iteration: 4000 / 8000 [ 50%] (Warmup)
Iteration: 4800 / 8000 [ 60%] (Sampling)
Iteration: 5600 / 8000 [ 70%] (Sampling)
Iteration: 6400 / 8000 [ 80%] (Sampling)
Iteration: 7200 / 8000 [ 90%] (Sampling)
Iteration: 8000 / 8000 [100%] (Sampling)
## Elapsed Time: 0.044 seconds (Warm-up)
## 0.056 seconds (Sampling)
## 0.1 seconds (Total)
##
## SAMPLING FOR MODEL 'mcmc4' NOW (CHAIN 4).
##
Iteration: 1 / 8000 [ 0%] (Warmup)
Iteration: 800 / 8000 [ 10%] (Warmup)
Iteration: 1600 / 8000 [ 20%] (Warmup)
Iteration: 2400 / 8000 [ 30%] (Warmup)
Iteration: 3200 / 8000 [ 40%] (Warmup)
Iteration: 4000 / 8000 [ 50%] (Warmup)
Iteration: 4800 / 8000 [ 60%] (Sampling)
Iteration: 5600 / 8000 [ 70%] (Sampling)
Iteration: 6400 / 8000 [ 80%] (Sampling)
Iteration: 7200 / 8000 [ 90%] (Sampling)
Iteration: 8000 / 8000 [100%] (Sampling)
## Elapsed Time: 0.049 seconds (Warm-up)
## 0.058 seconds (Sampling)
## 0.107 seconds (Total)
# output from stan
print(parSample, digits = 3)
## Inference for Stan model: mcmc4.
## 4 chains, each with iter=8000; warmup=4000; thin=4;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## m 1.142 0.005 0.281 0.603 0.966 1.143 1.314 1.718 3460 1
## s 0.853 0.004 0.223 0.539 0.695 0.812 0.960 1.401 3248 1
## lp__ -3.089 0.020 1.078 -5.967 -3.525 -2.745 -2.309 -2.027 3003 1
##
## Samples were drawn using NUTS(diag_e) at Thu Jun 12 19:40:52 2014.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
plot(parSample)
traceplot(parSample)
# get samples of parameters (representing integral over parameters)
par <- as.data.table(as.data.frame(parSample))
# rename col-names to match rest of code
setnames(par, c("m", "s"), c("mean", "sd"))
# add identifier & key
par[, `:=`(ipar, 1:nrow(par))]
## mean sd lp__ ipar
## 1: 0.9117 0.7539 -2.427 1
## 2: 1.2160 0.7887 -2.091 2
## 3: 1.1486 0.6616 -2.136 3
## 4: 1.0935 0.8216 -2.108 4
## 5: 1.0893 1.0788 -3.122 5
## ---
## 3996: 0.9953 0.8774 -2.376 3996
## 3997: 1.4264 0.8932 -2.850 3997
## 3998: 1.0841 0.8046 -2.079 3998
## 3999: 1.3197 1.2970 -4.335 3999
## 4000: 1.4307 0.8403 -2.774 4000
setkey(par, ipar)
date()
## [1] "Thu Jun 12 19:40:53 2014"
# simulate
dat <- as.data.table(my.sim(n, par))
# add identifier & key
dat[, `:=`(idat, 1:nrow(dat))]
## V1 V2 V3 V4 V5 V6 V7 V8
## 1: 0.6759 0.1281 1.1763 -0.1162 0.3082 1.4073 0.2407 0.8465
## 2: 0.3814 0.9294 1.8968 0.5768 1.6156 0.5547 0.5088 0.6246
## 3: 1.4246 0.2377 0.8317 1.6982 1.6532 -0.7317 -0.2761 0.6603
## 4: 0.8344 0.6536 0.1171 0.8270 -0.7400 0.1393 1.8257 1.0203
## 5: 1.9799 1.3858 0.7996 1.2500 1.0286 2.6721 -1.3002 -0.3902
## ---
## 3996: 1.4642 0.5346 0.9058 2.7452 1.0411 1.1023 2.4921 1.5314
## 3997: 0.9267 1.1563 1.6160 -0.2084 1.8354 2.0936 2.1870 0.4721
## 3998: 2.3781 1.4638 2.3969 2.1657 1.8666 0.3206 1.7595 0.7107
## 3999: 2.0417 0.9976 1.6501 -0.6007 -0.5129 1.7939 2.1831 0.2096
## 4000: 1.5867 2.4889 -0.3113 0.8280 2.0039 1.3176 0.5151 2.1727
## V9 V10 idat
## 1: 0.08350 -0.3149 1
## 2: 1.49599 0.5964 2
## 3: 2.01554 0.5644 3
## 4: 0.88062 1.1533 4
## 5: 1.09360 3.4431 5
## ---
## 3996: 2.18971 0.7507 3996
## 3997: 0.03271 1.9904 3997
## 3998: -0.40152 1.3898 3998
## 3999: -1.15258 0.8162 3999
## 4000: 0.37995 2.2445 4000
setkey(dat, idat)
# indices of data columns
iix <- match(grep("V", colnames(dat), value = TRUE), colnames(dat))
Given individual parameter vector as null hypothesis and complementary space as alternative hypothesis
# test a: test on difference of mean
fta <- function(par, dat) {
return(par[, abs(mean - mean(dat))])
}
# test b: student's t-test
ftb <- function(par, dat) {
return(par[, abs(mean - mean(dat))/sd(dat)])
}
# test c: likelihood ratio test ll ... log likelihood of null hypothesis,
# i.e., data given selected parameter mll ... log likelihood of alternative
# hypothesis, i.e., data and max. likl. parameter
ftc <- function(ll, mll) {
return(-2 * ll + 2 * mll)
}
# test d: posterior test - probability of parameter given data
ftd <- function(wn.par) {
return(1/wn.par)
}
# evaluate tests for simulated data
iix <- match(grep("V", colnames(dat), value = TRUE), colnames(dat))
llm[, `:=`(ta, fta(par, t(dat[idat, iix, with = FALSE]))), by = idat]
iix <- match(grep("V", colnames(dat), value = TRUE), colnames(dat))
llm[, `:=`(tb, ftb(par, t(dat[idat, iix, with = FALSE]))), by = idat]
llm[, `:=`(tc, ftc(ll, dat[idat, mll])), by = idat]
llm[, `:=`(td, ftd(wn.par)), by = idat]
# evaluate tests for observed data
par[, `:=`(tar, fta(par, x))]
par[, `:=`(tbr, ftb(par, x))]
par[, `:=`(tcr, ftc(ll, mll.obs$mll))]
par[, `:=`(tdr, ftd(1/nrow(par)))]
# evaluate tests using generic procedure
par <- par[llm[, list(p.a = sum(wn.sample * (ta > par[ipar, tar]))), by = ipar]]
par <- par[llm[, list(p.b = sum(wn.sample * (tb > par[ipar, tbr]))), by = ipar]]
par <- par[llm[, list(p.c = sum(wn.sample * (tc > par[ipar, tcr]))), by = ipar]]
par <- par[llm[, list(p.d = sum(wn.sample * (td > par[ipar, tdr]))), by = ipar]]
# student's t-test
par[, `:=`(p.t, pt(abs(mean - mean(x))/(sd(x)/sqrt(n)), df = n - 1, lower.tail = FALSE) *
2), by = ipar]
## ipar mean sd lp__ ll prior tar tbr tcr
## 1: 1 0.9117 0.7539 -2.427 -11.62 0.28246 0.22048 0.28245 0.8615
## 2: 2 1.2160 0.7887 -2.091 -11.28 0.23741 0.08377 0.10731 0.1889
## 3: 3 1.1486 0.6616 -2.136 -11.33 0.41314 0.01636 0.02096 0.2806
## 4: 4 1.0935 0.8216 -2.108 -11.30 0.19652 0.03866 0.04953 0.2237
## 5: 5 1.0893 1.0788 -3.122 -12.31 -0.07583 0.04292 0.05498 2.2524
## ---
## 3996: 3996 0.9953 0.8774 -2.376 -11.57 0.13077 0.13695 0.17544 0.7591
## 3997: 3997 1.4264 0.8932 -2.850 -12.04 0.11296 0.29418 0.37687 1.7072
## 3998: 3998 1.0841 0.8046 -2.079 -11.27 0.21744 0.04814 0.06167 0.1661
## 3999: 3999 1.3197 1.2970 -4.335 -13.52 -0.26006 0.18753 0.24024 4.6779
## 4000: 4000 1.4307 0.8403 -2.774 -11.96 0.17399 0.29850 0.38241 1.5561
## tdr p.a p.b p.c p.d p.t
## 1: 4000 0.3755 0.4285 0.6959 0.6556 0.3950
## 2: 4000 0.7333 0.7334 0.9138 0.9542 0.7421
## 3: 4000 0.9482 0.9557 0.8763 0.7285 0.9486
## 4: 4000 0.8726 0.8669 0.8959 0.9593 0.8790
## 5: 4000 0.8850 0.8549 0.3329 0.4244 0.8658
## ---
## 3996: 4000 0.6281 0.6033 0.7193 0.7917 0.5926
## 3997: 4000 0.2995 0.2629 0.4636 0.5136 0.2638
## 3998: 4000 0.8387 0.8428 0.9247 0.9679 0.8497
## 3999: 4000 0.6388 0.4430 0.1094 0.1311 0.4669
## 4000: 4000 0.2651 0.2569 0.4987 0.5206 0.2574
date()
## [1] "Thu Jun 12 19:42:02 2014"
ggplot(par, aes(x = mean, y = sd, colour = cut_number(p.a, n = 10))) + geom_point() +
geom_vline(aes(xintercept = mean(x))) + geom_hline(aes(yintercept = sd(x))) +
guides(colour = guide_legend("p-value")) + ggtitle("Naive Mean Test")
ggplot(par, aes(x = mean, y = sd, colour = cut_number(p.b, n = 10))) + geom_point() +
geom_vline(aes(xintercept = mean(x))) + geom_hline(aes(yintercept = sd(x))) +
guides(colour = guide_legend("p-value")) + ggtitle("T-Test")
ggplot(par, aes(x = mean, y = sd, colour = cut_number(p.c, n = 10))) + geom_point() +
geom_vline(aes(xintercept = mean(x))) + geom_hline(aes(yintercept = sd(x) *
sqrt((n - 1)/n))) + guides(colour = guide_legend("p-value")) + ggtitle("Likelihood-Ratio Test")
ggplot(par, aes(x = mean, y = sd, colour = cut_number(p.d, n = 10))) + geom_point() +
geom_vline(aes(xintercept = mean(x))) + geom_hline(aes(yintercept = sd(x))) +
guides(colour = guide_legend("p-value")) + ggtitle("Posterior Test")
ggplot(par[p.c > 0.95], aes(x = mean, y = sd, colour = cut_number(p.c, n = 10))) +
geom_point() + geom_vline(aes(xintercept = mean(x))) + geom_hline(aes(yintercept = sd(x) *
sqrt((n - 1)/n))) + guides(colour = guide_legend("p-value")) + ggtitle("Likelihood-Ratio Test")
ggplot(par, aes(x = p.c, y = p.d)) + geom_point()
ggplot(par, aes(x = mean)) + geom_point(aes(y = p.b), colour = "red") + geom_point(aes(y = p.t),
size = 0.1, colour = "yellow") + geom_vline(aes(xintercept = mean(x))) +
ylab("p-value") + ggtitle("T-Test")
ggplot(par, aes(x = tcr, y = p.c)) + geom_point(aes(y = 1 - pchisq(tcr, 2)),
colour = "yellow") + geom_point(alpha = 0.5) + xlab("Log Likelihood Ratio") +
ylab("p-value")
