网上看到一个学习资料，推荐下。参考资料：Chapter 7 Simulation-based Inference 来自于Book STAT160 R/RStudio Companion / 2021。
其主要介绍了4种推断方法
- One-Proportion Inference
- One-Mean Inference
- Two-proportion inference
- Two-mean inference
One-Proportion Inference
当你想测试某个人群比例不等于某个固定比例，那么你可以使用One-Proportion Inference。
比如你推断在一些教堂中女性的比例超过55%。这里的55%就是我们需要推断的proportion，hypotheses则是H0: pi=0.55和Ha: pi>0.55。然后我们发现有一个教堂中100名人员里面有62个女性，那么phat=62/100。此时我们是否有足够的证据推翻原假设？毕竟我们只是在一个教堂的抽样数据，所以我们得借助simulation，并假设原假设为真来计算假设检验的P值。因此我们借助do和rflip函数构建pi=0.55的1000个模拟trials
library(mosaic)
pi <- 0.55 # probability of success for each toss
n <- 100 # Number of times we toss the penny (sample size)
trials <- 1000 # Number of trials (number of samples)
observed <- 62 # Observed number of heads
phat = observed / n # p-hat - the observed proportion of heads
data.sim <- do(trials) * rflip(n, prob = pi)然后计算在这1000次模拟中，出现proportion大于phat的个数，除以总trials数就是我们想得到的"P值"。比较易于理解，由于是模拟所以得出来的，所以跟原本中的数值肯定是不会一样，但当模拟次数更大后，每次模拟的结果将会趋于一个较为稳定的值
pvalue <- sum(data.sim$prop >= phat) / trialsOne-Mean Inference
当你推测的值不是proportion而是一个mean值时，则可以考虑用One-Mean Inference。
比如你推测一辆汽车每加仑汽油的平均行驶里程不等于22英里。这里的22就是我们需要推断的mean，hypotheses则是H0: μ=22和Ha: μ≠22。然后我们观测到在数据集mtcars中每加仑汽油的平均行驶里程（mpg）的mean为20.09。这时我们可以对mtcars数据集进行重抽样来推断上述假设。
mu <- 22
observed <- mean(~mpg, data=mtcars, na.rm=T)
paste("Observed value for sample mean: ", observed)
trials <- 1000
samples <- do(trials) * mosaic::mean(~mpg, data=resample(mtcars))基于samples模拟数据，我们可以粗略计算下95%置信区间。
# Let's compute a 95% Confidence Interval
(ci <- quantile(samples$mean,c(0.025,0.975)))
# 2.5% 97.5%
# 18.18406 22.23453从置信区间可看出，其包含了我们H0假设的μ=22，所提可以初步判断原假设成立，即拒绝了汽车每加仑汽油的平均行驶里程不等于22这个推测。
接着根据重抽样的数据再次计算mpg的均值大于22的比例，由于是双侧假设，所以P值最终需要乘以2。
pvalue <- sum(samples$mean >= mu) / trials
paste("Two-sided p-value is", 2 * pvalue)
# [1] "Two-sided p-value is 0.088"Two-proportion inference
当你推测的是两个proportion之间是否有显著不同时，则可以考虑用Two-proportion inference。
比如你推测支持某项政策的女性比例与支持该政策的男性比例有不同。这里两个比例就是我们所需要比较的，hypotheses则是H0: π1=π2和Ha: π1≠π2。这时我们观测到某个样本里p1=62/100女性支持该政策，而男性则是p2=51/100，数据如下：
df <- rbind(
do(38) * data.frame(Group = "Men", Support = "no"),
do(62) * data.frame(Group = "Men", Support = "yes"),
do(49) * data.frame(Group = "Women", Support ="no"),
do(51) * data.frame(Group = "Women", Support = "yes")
)
(df.summary <- tally(Support ~ Group, data=df))接着计算两组之间proportion的差值，简单点就是0.62-0.51=0.11，或者
observed <- diffprop(Support ~ Group, data = df)
paste("Observed difference in proportions: " , round(observed,3))
# [1] "Observed difference in proportions: 0.11"然后使用打乱分组信息的方式做模拟，看看打乱后两组的差异的null distribution
trials <- 1000
null.dist <- do(trials) * diffprop(Support ~ shuffle(Group), data=df)
histogram( ~ diffprop, data= null.dist,
xlab = "Differences in proportions",
main = "Null distribution for differences in proportions",
v= observed)最后从null distribution中计算P值来确认是否能推翻原假设，也就是说假如H0假设是成立，那么差值大于0.11（或更大更极端的值）的概率是多少，是否很小（如小于0.05）。个人理解若P值很大，则可以推翻原假设。
p.value <- prop(~diffprop >= observed, data= null.dist)
paste(" One-sided p-value: ", round(p.value,3))
# [1] " One-sided p-value: 0.088"Two-mean inference
当你推测的是两个mean之间是否有显著不同时，则可以考虑用Two-mean inference
比如推测在长鳍金枪鱼（albacore）和黄鳍金枪鱼（yellowfin）中发现的汞的平均含量有不同。这里两个均值就是我们所需要比较的，hypotheses则是H0: μ1=μ2和Ha: μ1≠μ2。然后我们观测到在数据集tuna.txt中albacore mean为0.35763，而yellowfin mean为0.35444，两者的差值为-0.003。这时我们可以对tuna数据集进行打乱分组信息来模拟并推断上述假设。
df <- read.delim("http://citadel.sjfc.edu/faculty/ageraci/data/tuna.txt")
str(df)
favstats(~Mercury | Tuna, data=df)
observed <- diffmean(~Mercury | Tuna, data=df, na.rm = T)
paste("Observed difference in the means: ", round(observed, 3 ))
# [1] "Observed difference in the means: -0.003"接着就类似于Two-proportion inference，用shuffle函数打乱分组，模拟1000次，然后计算P值来判断当原假设成立前提下该值是否极端，最后看是否能推翻原假设。
trials <- 1000
null.dist <- do(trials) * diffmean(Mercury ~ shuffle(Tuna), data=df, na.rm = T)
pvalue <- prop(~ diffmean <= observed, data=null.dist)
paste("The one-sided p-value is ", round(pvalue,3))
# [1] "The one-sided p-value is 0.428"以上是对于参考资料的一个简单记录，simulation是一个非常有意思的方法，在临床试验中也较为常见，值得后续继续学习，本文所介绍的模拟是一个非常简单也易于理解的范例。
Reference
Chapter 7 Simulation-based Inference
Simulation-based inference with mosaic
MOSAIC R packages