Ancova仿真:功耗分析
我是一名大学生，正在上生态设计的统计课程，需要帮助找出如何在R中进行功率分析。我使用的是ANCOVA设计，没有交互作用；在这个假设的实验中，有两个花卉品种种植在独立的地块中，对于每个地块，收集的解释变量是土壤湿度，反应是花朵产量。
我模拟了一个名为sim.flowers的数据集，其中有花之间差异的α效应，斜率，x值的n，以及将在y向量中使用的正态分布的sigma (遵循y=α+ 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-值。我执行了以下操作，希望得到一个数据帧，其中有一列表示拒绝null的概率，这样我就可以将其绘制出来。
> 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")虽然我似乎理解功率分析所需代码的编写，但我在理解如何创建我应该创建的两个数字时遇到了麻烦。我不知道如何更改此代码，以使我生成的功率数字反映土壤湿度和花朵品种的影响。我非常感谢任何帮助我解决这个问题的人，我为我的长篇大论道歉，但我觉得这是我所问问题的必要背景。
回答 1
Stack Overflow用户
发布于 2018-07-04 21:48:59
您可以将您的结果与pwr.f2.test{pwr}进行比较
https://stackoverflow.com/questions/27372662
复制相似问题
腾讯云开发者
