私は、一種の遺伝的関連研究の検出力推定値を計算する Shiny アプリを持っています。ui.R は非常にシンプルで、server.R にはデータ フレームを提供する関数があります (いくつかのパラメーターがあるため、この関数を使用できないと思いreactive
ます)。
Gist へのリンクはhereです。実行するには:
library(shiny)
shiny:: runGist('5895082')
アプリは見積もりを正しく計算しますが、それに関して 2 つの質問があります。
output$powTable
最初の ? で、範囲内に含まれるすべての値を実際に表すことは可能sliderInput(n.cases)
ですか? 範囲の2つの極端な値を表しているようにしか見えません...何が間違っていますか?アプリの実行中にエラーが発生しました:
Error: Reading objects from shinyoutput object not allowed.
関数からデータ(反応性?)を渡しf()
てggplotにフィードするにはどうすればよいですか? 多くの試行錯誤の後、私は非常に迷っています。コードのどこでエラーが発生する可能性がありますか? 事前に多くの感謝を!
関数の元のコードはうまく機能します:(編集済み)
f <- function(ncases, p0, OR.cas.ctrl, Nh, sig.level) {
num.cases <- ncases
p0 <- p0
Nh <- Nh
OR.cas.ctrl <- OR.cas.ctrl
sig.level <- sig.level
# Parameters related to sig.level, from [Table 2] of Samuels et al.
# For 90% power and alpha = .05, Nscaled = 8.5
if (sig.level == 0.05){
A <- -28 # Parameter A for alpha=.05
x0 <- 2.6 # Parameter x0 for alpha=.05
d <- 2.4 # Parameter d for alpha=.05
}
if (sig.level == 0.01){
A <- -13 # Parameter A for alpha=.01
x0 <- 5 # Parameter x0 for alpha=.01
d <- 2.5 # Parameter d for alpha=.01
}
if (sig.level == 0.001){
A <- -7 # Parameter A for alpha=.001
x0 <- 7.4 # Parameter x0 for alpha=.001
d <- 2.8 # Parameter d for alpha=.001
}
out.pow <- NULL # initialize vector
for(ncases in ncases){
OR.ctrl.cas <- 1 / OR.cas.ctrl # 1. CALCULATE P1 FROM A PREDEFINED P0, AND A DESIRED OR
OR <- OR.ctrl.cas
bracket.pw <- p0 / (OR - OR*p0) # obtained after isolating p1 in OR equation [3].
p1 <- bracket.pw / (1 + bracket.pw)
Nh037 <- Nh^0.37 # 2. CALCULATE NSCALED
num.n <- num.cases*((p1-p0)^2)
den.n <- (p1*(1-p1) + p0*(1-p0))*Nh037
Nscaled <- num.n/den.n
num.power <- A - 100 # 3. CALCULATE POWER
den.power <- 1 + exp((Nscaled - x0)/d)
power <- 100 + (num.power/den.power) # The power I have to detect a given OR with my data, at a given alpha
}
OR <- OR.cas.ctrl
out.pow <- data.frame(num.cases, Nh, Nscaled, p0, OR, sig.level, power)
out.pow
}
mydata <- f(ncases=seq(50,1000, by=50), 0.4, 2.25, 11, 0.05)
mydata
library(ggplot2)
print(ggplot(data = mydata, aes(num.cases, power)) +
theme_bw() +
theme(text=element_text(family="Helvetica", size=12)) +
labs(title = "Ad-hoc power for haplogroup") +
scale_color_brewer(palette = "Dark2", guide = guide_legend(reverse=TRUE)) +
xlab("number of cases/controls") +
ylab("power") +
scale_x_log10() +
geom_line(alpha=0.8, size=0.2) +
geom_point(aes(shape = factor(OR)), colour="black"))