Ancova 模拟：功率分析

Question

我是一名正在上生态设计统计课程的大学生，需要帮助弄清楚如何在 R 中进行功效分析。我正在使用没有交互的 ANCOVA 设计；在这个假设的实验中，在独立的地块中种植了两种花卉品种，对于每个地块，收集的解释变量是土壤水分，响应是花卉产量。

我模拟了一个名为 sim.flowers 的数据集，其中包含花朵之间差异的 alpha 效应、斜率、我的 x 值的 n 以及我将在我的 y 向量中使用的正态分布的 sigma（以下y = alpha + beta0 + beta1X 的 ancova 模型；为简单起见，我将截距设为 0) 见下文：

sim.flowers <- function(alpha,slope,n,sigma) {
x <-runif(2*n, min = -1, max = 1) 
flower.effects <- rep(c(0,alpha),each=n) #there are two different flower varieties and I gave them a true difference of 1. 
y <- flower.effects + slope * x  + rnorm(2*n, 0, sigma)
data.frame(x=x,y=y,flower.effects = flower.effects)
}

我对此进行了测试，它成功了，它给了我一个带有 X a Y 和花卉效果列的数据集

> test1 <- sim.oats(1,0.5,3,0.3)
> test1
        x           y        flower.effects
1 -0.99913780 -0.31373866           0
2 -0.38610391  0.41070965           0
3 -0.58308522  0.07426254           0
4 -0.35900237  0.36395132           1
5 -0.07296464  1.29149447           1
6  0.18575996  0.85001847           1

我们的目标是创建一个显示检测花卉品种效应的能力的图形，对于每个花卉品种处理的不同重复次数，我被告知应该只选择土壤水分效应的 1 个值。

以及显示检测土壤水分效应的功效的图，对于每个花卉品种处理，不同的重复次数具有不同的线条，选择花卉品种的 1 个值

为了开始到达这里，我在 for 循环中运行了线性回归，以便我可以提取 p 值并能够绘制一个幂图，将拒绝空值的概率设置为 0.5 我的代码低于

> sim.flowers.many <- function(alpha,slope,n,sigma,numsimulations){
+   pvals <-numeric(numsimulations)
+   for(i in 1:numsimulations){
+       thisdat <-sim.flowers(alpha,slope,n,sigma)
+       thisfit <-lm(y~x,thisdat)
+     pvals[i]<-coefficients(summary(thisfit))['x','Pr(>|t|)']
+   }
+   return(pvals)
+ }
> sim.flowers.many( alpha = 1,slope = 0.5, n = 3, sigma = 6, numsimulations =3)
[1] 0.7662218 0.4454654 0.2414637

我似乎得到了 p 值。我做了以下操作，期望我会得到一个数据框，其中有一列说明拒绝空值的概率，以便我可以将其绘制出来。

> determine.power <- function(true.slopes){
+     true.slopes <-0:3
+     out<- data.frame(true.slopes,prob.reject.null=NA)
+       for(i in 1:length(true.slopes)){
+       thesepvals <- sim.flowers.many(alpha = 1,slope = 0.5, n = 3, sigma = 6, numsimulations =3)
+       out[i,2] <- mean(thesepvals < 0.05)
+         }
+     return(out)
+     }

我得到了这个作为输出：

> determine.power(true.slopes=0.5)
    true.slopes prob.reject.null
1           0        0.0000000
2           1        0.0000000
3           2        0.3333333
4           3        0.3333333

我想我可以这样画出来：

 power.results <- determine.power(seq(0, 2, by = 0.2))
 plot(prob.reject.null ~ true.slopes, data = power.results, main = "Power analysis", xlab     = "true slope", 
ylab = "Prob. reject null", ylim = c(0,1), type = "b", col = "slateblue4") 
 grid(col = "hotpink")

虽然我似乎理解了编写功率分析所需的代码，但我无法理解如何创建我应该创建的两个图形。我不知道如何更改此代码，以便我生成的功率数字反映土壤水分和花卉品种的影响。对于解决这个问题的任何帮助，我将不胜感激，我为篇幅道歉，但我认为这是我所问问题的必要背景。

Answer 1

您可以将您的结果与pwr.f2.test{pwr}进行比较