DALS005-统计推断(Inference)04-Monte

---

title: DALS005-统计推断(Inference)04-Monte Carlo
date: 2019-08-05 12:0:00
type: "tags"
tags:

• 统计推断
  categories:
• 生物统计

---

前言

计算机可以生成伪随机数(pseudo-random numbers)，使用这些伪随机数可以用于模拟来自真实世界的随机变量。这可以让我们利用计算机来，而不是理论分析或推导来研究一些随机变量的特异。使用这种概念的一个非常有用的地方就在于，我们可以创建一些模拟数据来测试来验证我们的一些想法，而不必在实验室中做一些实验。

模拟数据也可以用于检测一些理论或分析结果的可靠性。此外，在统计学中，我们使用的很多结果是渐时性的：只当样本达到无穷大时它们才有可能成立。在实际中，我们无法收集所有的样本来验证理论与实际的关系如何。有时候我们就需要一些模拟数据。

Monte Carlo模拟

现在我们使用Mote Carlo模拟来比较一下在不同样本数目中，CLT与t分布的接近程度，如下所示：

library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% dplyr::select(Bodyweight) %>% unlist

现在构建一个函数，这个函数用于生成零假设下，样本数目n的t统计量，如下所示：

library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% dplyr::select(Bodyweight) %>% unlist

ttestgenerator <- function(n){
  # note that here we have a false "high fat" group where we actually
  # sample from the nonsmokers. this is because we are modeling the *null*
  cases <- sample(controlPopulation, n)
  controls <- sample(controlPopulation,n)
  tstat <- (mean(cases) - mean(controls))/
    sqrt(var(cases)/n + var(controls)/n)
  return(tstat)
}
ttests <- replicate(1000,ttestgenerator(10))
hist(ttests)

image

这个t分布是否非常接近于正态分布？我们可以使用qq图(quantile-quantile)来看一下，如下所示：

qqnorm(ttests)
abline(0,1)

image

从上图我们可以看出来，这个数据集非常接近于正态分布。如果现在我们将样本数改为3，qq图如下所示：

ttests <- replicate(1000, ttestgenerator(3))
qqnorm(ttests)
abline(0,1)

image

从上图可以看出来，数据的高分位数(large quantiles)区，也就是统计学家称的尾部(tails)，比较大，这个尾部区就是上图中直线左侧的下方数据，以及右侧直线上方的数据。在前面的部分我们提到，即使总体服从正态分布，那么小样本则是更加近似地服从t分布，下面我们来的模拟一下：

ps <- (seq(0,999)+0.5)/1000
qqplot(qt(ps,df=2*3-2),ttests,xlim=c(-6,6),ylim=c(-6,6))
abline(0,1)

image

从上面两张图可以看出来，数据更接近于t分布，而不是正态分布，不过虽然接近于t分布，但是这种接近也并非完全，这是由于原始数据也并不是完美的服从正态分布，如下所示：

qqnorm(controlPopulation)
qqline(controlPopulation)

image

观测值的参数化模拟(Parametric Simulations for the Observations)

上述我们模拟随机变量的方法叫Monte Carlo模拟。我们利用既有的总体数据，随机抽取了一些样本。在实际分析中，我们无法获取整个总体。但是，当我们想在实际分析过程中使用Monte Carlo模拟时，一种比较典型的做法是，假设一个参数分布(parametric distribution)，并根据此分布生成一个总体，这种手段叫参数模拟(parametric simulation)。这就意味着，我们从真实的数据（这里指的是样本的均值与标准差）中提取参数，将这些参数添加一个模型（这里指的是正态分布）。这是Monte Carlo模拟最普遍的手段。

回到小鼠的体重案例上来，我们可以利用我们现在有的知识，例如小鼠的体重通常是24g，SD为3.5g，小鼠的体重用了近服从正态分布，我们可以产生一个总体：

controls<- rnorm(5000, mean=24, sd=3.5)

当我们生成了这批数据后，我们可以重复再生成几次。此时，我们并不使用sample()函数来进行抽样，而是使用随机数来生成一些服从正态分布的数据，前面提到的ttestgenerator()函数更改如下：

ttestgenerator <- function(n, mean=24, sd=3.5){
  cases <- rnorm(n, mean,sd)
  controls <- rnorm(n, mean, sd)
  tstat <- (mean(cases)-mean(controls))/
    sqrt(var(cases)/n + var(controls)/n)
  return(tstat)
}

练习

P88