lvr-meta.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Meta-analysis on log-variance ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulation parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "suppressPackageStartupMessages(library(meta))\n",
    "\n",
    "Niter <- 10000          # Number of iterations for bootstrap simulation\n",
    "Nstud <- 10             # Number of studies in the meta-analysis\n",
    "ss1 <- 20:40            # Different possible number of cases\n",
    "ss2 <- 5:30             # Different possible number of controls\n",
    "mu1 <- 10               # True mean for cases\n",
    "mu2 <- 10               # True mean for controls\n",
    "sigma1 <- 1.2           # True standard deviation for cases\n",
    "sigma2 <- 1             # True standard deviation for controls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Samples generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n1 <- sample(ss1, Nstud, rep=T)     # Number of cases for each study\n",
    "n2 <- sample(ss2, Nstud, rep=T)     # Number of controls for each study\n",
    "\n",
    "ind1 <- lapply(n1, rnorm, mu1, sigma1)      # Cases generation\n",
    "ind2 <- lapply(n2, rnorm, mu2, sigma2)      # Controls generation\n",
    "sd1 <- sapply(ind1, sd)                     # Estimated standard deviations for cases of each study\n",
    "sd2 <- sapply(ind2, sd)                     # Estimated standard deviations for controls of each study"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### F-test for comparison two variances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>study</th><th scope=col>cases</th><th scope=col>controls</th><th scope=col>estimate</th><th scope=col>lower.95</th><th scope=col>upper.95</th><th scope=col>pval</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td> 1        </td><td>23        </td><td>24        </td><td>1.8848013 </td><td>0.8105226 </td><td>4.418297  </td><td>0.13867008</td></tr>\n",
       "\t<tr><td> 2        </td><td>23        </td><td>12        </td><td>1.3381884 </td><td>0.4185702 </td><td>3.541987  </td><td>0.63052687</td></tr>\n",
       "\t<tr><td> 3        </td><td>24        </td><td>24        </td><td>1.9907306 </td><td>0.8611765 </td><td>4.601854  </td><td>0.10575130</td></tr>\n",
       "\t<tr><td> 4        </td><td>28        </td><td>15        </td><td>0.8789126 </td><td>0.3187168 </td><td>2.104866  </td><td>0.74520592</td></tr>\n",
       "\t<tr><td> 5        </td><td>34        </td><td>22        </td><td>0.9191367 </td><td>0.4020617 </td><td>1.959113  </td><td>0.80907985</td></tr>\n",
       "\t<tr><td> 6        </td><td>28        </td><td> 9        </td><td>0.9869077 </td><td>0.2519044 </td><td>2.671950  </td><td>0.89910365</td></tr>\n",
       "\t<tr><td> 7        </td><td>23        </td><td>29        </td><td>1.5539479 </td><td>0.7061529 </td><td>3.557181  </td><td>0.26937288</td></tr>\n",
       "\t<tr><td> 8        </td><td>23        </td><td>21        </td><td>2.0175320 </td><td>0.8290070 </td><td>4.819849  </td><td>0.11963258</td></tr>\n",
       "\t<tr><td> 9        </td><td>20        </td><td>17        </td><td>2.4768702 </td><td>0.9180296 </td><td>6.416713  </td><td>0.07229135</td></tr>\n",
       "\t<tr><td>10        </td><td>26        </td><td>12        </td><td>2.8771494 </td><td>0.9100170 </td><td>7.366377  </td><td>0.07043425</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllllll}\n",
       " study & cases & controls & estimate & lower.95 & upper.95 & pval\\\\\n",
       "\\hline\n",
       "\t  1         & 23         & 24         & 1.8848013  & 0.8105226  & 4.418297   & 0.13867008\\\\\n",
       "\t  2         & 23         & 12         & 1.3381884  & 0.4185702  & 3.541987   & 0.63052687\\\\\n",
       "\t  3         & 24         & 24         & 1.9907306  & 0.8611765  & 4.601854   & 0.10575130\\\\\n",
       "\t  4         & 28         & 15         & 0.8789126  & 0.3187168  & 2.104866   & 0.74520592\\\\\n",
       "\t  5         & 34         & 22         & 0.9191367  & 0.4020617  & 1.959113   & 0.80907985\\\\\n",
       "\t  6         & 28         &  9         & 0.9869077  & 0.2519044  & 2.671950   & 0.89910365\\\\\n",
       "\t  7         & 23         & 29         & 1.5539479  & 0.7061529  & 3.557181   & 0.26937288\\\\\n",
       "\t  8         & 23         & 21         & 2.0175320  & 0.8290070  & 4.819849   & 0.11963258\\\\\n",
       "\t  9         & 20         & 17         & 2.4768702  & 0.9180296  & 6.416713   & 0.07229135\\\\\n",
       "\t 10         & 26         & 12         & 2.8771494  & 0.9100170  & 7.366377   & 0.07043425\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "study | cases | controls | estimate | lower.95 | upper.95 | pval | \n",
       "|---|---|---|---|---|---|---|---|---|---|\n",
       "|  1         | 23         | 24         | 1.8848013  | 0.8105226  | 4.418297   | 0.13867008 | \n",
       "|  2         | 23         | 12         | 1.3381884  | 0.4185702  | 3.541987   | 0.63052687 | \n",
       "|  3         | 24         | 24         | 1.9907306  | 0.8611765  | 4.601854   | 0.10575130 | \n",
       "|  4         | 28         | 15         | 0.8789126  | 0.3187168  | 2.104866   | 0.74520592 | \n",
       "|  5         | 34         | 22         | 0.9191367  | 0.4020617  | 1.959113   | 0.80907985 | \n",
       "|  6         | 28         |  9         | 0.9869077  | 0.2519044  | 2.671950   | 0.89910365 | \n",
       "|  7         | 23         | 29         | 1.5539479  | 0.7061529  | 3.557181   | 0.26937288 | \n",
       "|  8         | 23         | 21         | 2.0175320  | 0.8290070  | 4.819849   | 0.11963258 | \n",
       "|  9         | 20         | 17         | 2.4768702  | 0.9180296  | 6.416713   | 0.07229135 | \n",
       "| 10         | 26         | 12         | 2.8771494  | 0.9100170  | 7.366377   | 0.07043425 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "   study cases controls estimate  lower.95  upper.95 pval      \n",
       "1   1    23    24       1.8848013 0.8105226 4.418297 0.13867008\n",
       "2   2    23    12       1.3381884 0.4185702 3.541987 0.63052687\n",
       "3   3    24    24       1.9907306 0.8611765 4.601854 0.10575130\n",
       "4   4    28    15       0.8789126 0.3187168 2.104866 0.74520592\n",
       "5   5    34    22       0.9191367 0.4020617 1.959113 0.80907985\n",
       "6   6    28     9       0.9869077 0.2519044 2.671950 0.89910365\n",
       "7   7    23    29       1.5539479 0.7061529 3.557181 0.26937288\n",
       "8   8    23    21       2.0175320 0.8290070 4.819849 0.11963258\n",
       "9   9    20    17       2.4768702 0.9180296 6.416713 0.07229135\n",
       "10 10    26    12       2.8771494 0.9100170 7.366377 0.07043425"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Fr <- sd1^2 / sd2^2                                                 # Fisher ratios\n",
    "df1 <- n1 - 1                                                       # Degrees of freedom for F-test\n",
    "df2 <- n2 - 1\n",
    "f.pvalues <- 1 - abs(pf(Fr, df1, df2) - 0.5) * 2                    # F-test p-values\n",
    "alpha <- .05\n",
    "ci <- cbind(Fr/qf(1-alpha/2, df1, df2), Fr/qf(alpha/2, df1, df2))   # 95% Confidence intervals for variance ratios\n",
    "data.frame(study=1:Nstud, cases=n1, controls=n2, estimate=Fr, lower.95=ci[,1], upper.95=ci[,2], pval=f.pvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Meta-analysis on ratios of variances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variance ratio estimate (real value: 1.44):"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>method</th><th scope=col>estimate</th><th scope=col>lower.95</th><th scope=col>upper.95</th><th scope=col>pval</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>fixed      </td><td>1.239225   </td><td>1.07264    </td><td>1.431681   </td><td>0.003591341</td></tr>\n",
       "\t<tr><td>random     </td><td>1.239225   </td><td>1.07264    </td><td>1.431681   </td><td>0.003591341</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllll}\n",
       " method & estimate & lower.95 & upper.95 & pval\\\\\n",
       "\\hline\n",
       "\t fixed       & 1.239225    & 1.07264     & 1.431681    & 0.003591341\\\\\n",
       "\t random      & 1.239225    & 1.07264     & 1.431681    & 0.003591341\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "method | estimate | lower.95 | upper.95 | pval | \n",
       "|---|---|\n",
       "| fixed       | 1.239225    | 1.07264     | 1.431681    | 0.003591341 | \n",
       "| random      | 1.239225    | 1.07264     | 1.431681    | 0.003591341 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  method estimate lower.95 upper.95 pval       \n",
       "1 fixed  1.239225 1.07264  1.431681 0.003591341\n",
       "2 random 1.239225 1.07264  1.431681 0.003591341"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "TE <- log(Fr) - digamma(df1/2) + digamma(df2/2) + log(df1/df2)\n",
    "TEse <- sqrt(trigamma(df1/2) + trigamma(df2/2))\n",
    "meta <- metagen(TE, TEse)\n",
    "cat(sprintf(\"Variance ratio estimate (real value: %g):\", sigma1^2/sigma2^2))\n",
    "with(meta, data.frame(method=c(\"fixed\", \"random\"), estimate=exp(c(TE.fixed, TE.random)/2),\n",
    "        lower.95=exp(c(lower.fixed, lower.random)/2),\n",
    "        upper.95=exp(c(upper.fixed, upper.random)/2), pval=c(pval.fixed, pval.random)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "metas <- replicate(Niter, {\n",
    "\tn1 <- sample(ss1, Nstud, rep=T)\n",
    "\tn2 <- sample(ss2, Nstud, rep=T)\n",
    "\t\n",
    "\tind1 <- lapply(n1, rnorm, mu1, sigma1)\n",
    "\tind2 <- lapply(n2, rnorm, mu2, sigma2)\n",
    "\tsd1 <- sapply(ind1, sd)\n",
    "\tsd2 <- sapply(ind2, sd)\n",
    "\t\n",
    "\tFr <- sd1^2 / sd2^2\n",
    "\tdf1 <- n1 - 1\n",
    "\tdf2 <- n2 - 1\n",
    "    \n",
    "\tTE1 <- log(Fr) - digamma(df1/2) + digamma(df2/2) + log(df1/df2)     # Correcting for bias\n",
    "\tTE2 <- log(Fr)                                                      # Not correcting for bias\n",
    "\tTEse <- sqrt(trigamma(df1/2) + trigamma(df2/2))\n",
    "\tlist(nobias=metagen(TE1, TEse, n.e=n1, n.c=n2), bias=metagen(TE2, TEse, n.e=n1, n.c=n2))\n",
    "}, simplify=F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Distribution of estimations\n",
    "<b>mean:</b> mean log-variance ratio estimated by meta-analysis<br>\n",
    "<b>sd:</b> standard deviation of the estimations, which corresponds to the standard error of one estimation<br>\n",
    "<b>sem:</b> standard error for the mean of the estimations<br>\n",
    "<b>p.value:</b> p-value for the existence of a bias in the estimations<br>\n",
    "<b>fixed:</b> fixed effect model<br>\n",
    "<b>random:</b> random effects model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of cases by meta-analysis: 300.098 ± 19.5001\n",
      "Total number of cases by meta-analysis: 174.579 ± 23.6432\n",
      "\n",
      "Real log-variance ratio: 0.364643\n",
      "\n",
      "Combined estimations with bias correction:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>sem</th><th scope=col>p.value</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>fixed</th><td>0.3622764  </td><td>0.1470438  </td><td>0.001470438</td><td>0.1075013  </td></tr>\n",
       "\t<tr><th scope=row>random</th><td>0.3634180  </td><td>0.1480556  </td><td>0.001480556</td><td>0.4079570  </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llll}\n",
       "  & mean & sd & sem & p.value\\\\\n",
       "\\hline\n",
       "\tfixed & 0.3622764   & 0.1470438   & 0.001470438 & 0.1075013  \\\\\n",
       "\trandom & 0.3634180   & 0.1480556   & 0.001480556 & 0.4079570  \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | mean | sd | sem | p.value | \n",
       "|---|---|\n",
       "| fixed | 0.3622764   | 0.1470438   | 0.001470438 | 0.1075013   | \n",
       "| random | 0.3634180   | 0.1480556   | 0.001480556 | 0.4079570   | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "       mean      sd        sem         p.value  \n",
       "fixed  0.3622764 0.1470438 0.001470438 0.1075013\n",
       "random 0.3634180 0.1480556 0.001480556 0.4079570"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Combined estimations without bias correction:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>mean</th><th scope=col>sd</th><th scope=col>sem</th><th scope=col>p.value</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>fixed</th><td>0.3939271  </td><td>0.1473636  </td><td>0.001473636</td><td>0          </td></tr>\n",
       "\t<tr><th scope=row>random</th><td>0.3974276  </td><td>0.1490018  </td><td>0.001490018</td><td>0          </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llll}\n",
       "  & mean & sd & sem & p.value\\\\\n",
       "\\hline\n",
       "\tfixed & 0.3939271   & 0.1473636   & 0.001473636 & 0          \\\\\n",
       "\trandom & 0.3974276   & 0.1490018   & 0.001490018 & 0          \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | mean | sd | sem | p.value | \n",
       "|---|---|\n",
       "| fixed | 0.3939271   | 0.1473636   | 0.001473636 | 0           | \n",
       "| random | 0.3974276   | 0.1490018   | 0.001490018 | 0           | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "       mean      sd        sem         p.value\n",
       "fixed  0.3939271 0.1473636 0.001473636 0      \n",
       "random 0.3974276 0.1490018 0.001490018 0      "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "metas1 <- lapply(metas, \"[[\", \"nobias\")\n",
    "metas2 <- lapply(metas, \"[[\", \"bias\")\n",
    "Ne <- colSums(sapply(metas1, '[[', 'n.e'))\n",
    "Nc <- colSums(sapply(metas1, '[[', 'n.c'))\n",
    "lvr <- 2*log(sigma1/sigma2)\n",
    "cat(sprintf(\"Total number of cases by meta-analysis: %g ± %g\\n\", mean(Ne), sd(Ne)))\n",
    "cat(sprintf(\"Total number of cases by meta-analysis: %g ± %g\\n\\n\", mean(Nc), sd(Nc)))\n",
    "cat(sprintf(\"Real log-variance ratio: %g\\n\\n\", lvr))\n",
    "cat(\"Combined estimations with bias correction:\\n\")\n",
    "TEs1.fixed <- sapply(metas1, \"[[\", \"TE.fixed\")\n",
    "TEs1.random <- sapply(metas1, \"[[\", \"TE.random\")\n",
    "TEs1.mean <- c(mean(TEs1.fixed), mean(TEs1.random))\n",
    "TEs1.sd <- c(sd(TEs1.fixed), sd(TEs1.random))\n",
    "TEs1.sem <- TEs1.sd / sqrt(length(TEs1.fixed))\n",
    "TEs1.p <- 1-2*abs(pnorm(TEs1.mean, lvr, TEs1.sem)-.5)\n",
    "data.frame(mean = TEs1.mean, sd=TEs1.sd, sem=TEs1.sem, p.value=TEs1.p, row.names=c(\"fixed\", \"random\"))\n",
    "cat(\"Combined estimations without bias correction:\\n\")\n",
    "TEs2.fixed <- sapply(metas2, \"[[\", \"TE.fixed\")\n",
    "TEs2.random <- sapply(metas2, \"[[\", \"TE.random\")\n",
    "TEs2.mean <- c(mean(TEs2.fixed), mean(TEs2.random))\n",
    "TEs2.sd <- c(sd(TEs2.fixed), sd(TEs2.random))\n",
    "TEs2.sem <- TEs2.sd / sqrt(length(TEs2.fixed))\n",
    "TEs2.p <- 1-2*abs(pnorm(TEs2.mean, lvr, TEs2.sem)-.5)\n",
    "data.frame(mean = TEs2.mean, sd=TEs2.sd, sem=TEs2.sem, p.value=TEs2.p, row.names=c(\"fixed\", \"random\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistical power of estimations (rate of log-variance ratios estimated significantly different from 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>nobias</th><th scope=col>bias</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>fixed</th><td>0.6940</td><td>0.7619</td></tr>\n",
       "\t<tr><th scope=row>random</th><td>0.6285</td><td>0.7084</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & nobias & bias\\\\\n",
       "\\hline\n",
       "\tfixed & 0.6940 & 0.7619\\\\\n",
       "\trandom & 0.6285 & 0.7084\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | nobias | bias | \n",
       "|---|---|\n",
       "| fixed | 0.6940 | 0.7619 | \n",
       "| random | 0.6285 | 0.7084 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "       nobias bias  \n",
       "fixed  0.6940 0.7619\n",
       "random 0.6285 0.7084"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "power1.fixed <- mean(sapply(metas1, \"[[\", \"pval.fixed\") <= alpha)\n",
    "power1.random <- mean(sapply(metas1, \"[[\", \"pval.random\") <= alpha)\n",
    "power2.fixed <- mean(sapply(metas2, \"[[\", \"pval.fixed\") <= alpha)\n",
    "power2.random <- mean(sapply(metas2, \"[[\", \"pval.random\") <= alpha)\n",
    "data.frame(nobias=c(power1.fixed, power1.random), bias=c(power2.fixed, power2.random), row.names=c(\"fixed\", \"random\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Frequency of incorrect estimations (Rate of estimated 95% confidence intervals not including real value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>nobias</th><th scope=col>bias</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>fixed</th><td>0.0507</td><td>0.0536</td></tr>\n",
       "\t<tr><th scope=row>random</th><td>0.0384</td><td>0.0406</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & nobias & bias\\\\\n",
       "\\hline\n",
       "\tfixed & 0.0507 & 0.0536\\\\\n",
       "\trandom & 0.0384 & 0.0406\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| <!--/--> | nobias | bias | \n",
       "|---|---|\n",
       "| fixed | 0.0507 | 0.0536 | \n",
       "| random | 0.0384 | 0.0406 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "       nobias bias  \n",
       "fixed  0.0507 0.0536\n",
       "random 0.0384 0.0406"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lvr <- 2*log(sigma1/sigma2)\n",
    "fie1.fixed <- mean(sapply(metas1, \"[[\", \"lower.fixed\") > lvr | sapply(metas1, \"[[\", \"upper.fixed\") < lvr)\n",
    "fie1.random <- mean(sapply(metas1, \"[[\", \"lower.random\") > lvr | sapply(metas1, \"[[\", \"upper.random\") < lvr)\n",
    "fie2.fixed <- mean(sapply(metas2, \"[[\", \"lower.fixed\") > lvr | sapply(metas2, \"[[\", \"upper.fixed\") < lvr)\n",
    "fie2.random <- mean(sapply(metas2, \"[[\", \"lower.random\") > lvr | sapply(metas2, \"[[\", \"upper.random\") < lvr)\n",
    "data.frame(nobias=c(fie1.fixed, fie1.random), bias=c(fie2.fixed, fie2.random), row.names=c(\"fixed\", \"random\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.3.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}