Bonferroni and conditional coverage

Bonferroni doesn't provide conditional coverage.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "NULL"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "simulate = function(n=200,\n",
    "                    p=100,\n",
    "                    s=10,\n",
    "                    signal_range=c(0.5, 1),\n",
    "                    rho=0.5,\n",
    "                    q=0.2,\n",
    "                    level=0.9) {\n",
    "    \n",
    "    beta = rep(0, p)\n",
    "    beta[sample(1:p, s, replace=FALSE)] = seq(signal_range[1], signal_range[2], length=s) / sqrt(n)\n",
    "    \n",
    "    X0 = matrix(rnorm(n * (p+1)), n, p+1)\n",
    "    X = X0[,2:(p+1)] + 0.5 * X0[,1:p]\n",
    "    X = scale(X, TRUE, TRUE)\n",
    "    y = X %*% beta + rnorm(n)\n",
    "    \n",
    "    M = lm(y ~ X - 1)\n",
    "    pval = coef(summary(M))[,4]\n",
    "    selected = p.adjust(pval, method='BH') < q\n",
    "    \n",
    "    alpha = 1 - level\n",
    "    bonf_alpha = alpha / p\n",
    "    bonf_level = 1 - bonf_alpha\n",
    "    \n",
    "    # BY does bonferroni with num selected\n",
    "    \n",
    "    if (sum(selected) > 0) {\n",
    "       BY_alpha = alpha / sum(selected)\n",
    "       BY_level = 1 - BY_alpha\n",
    "       selected_intervals = confint(M, level=bonf_level)[selected,,drop=FALSE]\n",
    "       selected_BY = confint(M, level=BY_level)[selected,,drop=FALSE]\n",
    "       targets = beta[selected]\n",
    "       bonf_coverage = (selected_intervals[,1] < targets) * (selected_intervals[,2] > targets)\n",
    "       BY_coverage = (selected_BY[,1] < targets) * (selected_BY[,2] > targets)\n",
    "       coverage = data.frame(bonferroni=bonf_coverage, BY=BY_coverage)\n",
    "       FCP = data.frame(bonferroni=1 - mean(bonf_coverage), BY=1 - mean(BY_coverage))\n",
    "       results = list(coverage=coverage, FCP=FCP)\n",
    "    } else { \n",
    "        results = NULL\n",
    "    }\n",
    "    return(results)\n",
    "}\n",
    "\n",
    "simulate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Raw proportion of intervals covering their target\"\n",
      "bonferroni         BY \n",
      " 0.8243993  0.1700555 \n",
      "[1] \"False coverage rate\"\n",
      "[1] 0.2570205 0.9628271\n",
      "bonferroni         BY \n",
      "0.03318854 0.15685952 \n"
     ]
    }
   ],
   "source": [
    "coverage = c()\n",
    "FCP = c()\n",
    "zero_reported = c()\n",
    "\n",
    "for (i in 1:1000) {\n",
    "    results = simulate()\n",
    "    if (!is.null(results)) {\n",
    "       coverage = rbind(coverage, results$coverage)\n",
    "       FCP = rbind(FCP, results$FCP)\n",
    "    }\n",
    "    zero_reported = c(zero_reported, is.null(results))\n",
    "}\n",
    "print('Raw proportion of intervals covering their target')\n",
    "print(apply(coverage, 2, mean))\n",
    "\n",
    "print('False coverage rate')\n",
    "print(c(mean(FCP[,'bonferroni']), mean(FCP[,'BY'])))\n",
    "\n",
    "# Bonferroni should control at 1 - level the following\n",
    "print((1 - apply(coverage, 2, mean)) * (1 - mean(zero_reported)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "all,-slideshow"
  },
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}