Leistungsanalyse für ordinale logistische Regression

Ich bin auf der Suche nach einem Programm (in R oder SAS oder eigenständig, wenn es kostenlos oder kostengünstig ist), das eine Leistungsanalyse für die ordinale logistische Regression durchführt.

Ich bevorzuge Leistungsanalysen jenseits der Grundlagen durch Simulation. Bei vorgefertigten Paketen bin ich mir nie ganz sicher, welche Annahmen getroffen werden.

Das Simulieren für Leistung ist mit R ganz einfach (und erschwinglich).

1. Entscheiden Sie, wie Ihre Daten Ihrer Meinung nach aussehen sollen und wie Sie sie analysieren werden
2. Schreiben Sie eine Funktion oder eine Gruppe von Ausdrücken, die die Daten für eine bestimmte Beziehung und Stichprobengröße simulieren und die Analyse durchführen (eine Funktion ist vorzuziehen, da Sie die Stichprobengröße und die Parameter in Argumente umwandeln können, um das Ausprobieren verschiedener Werte zu erleichtern). Die Funktion oder der Code sollte den p-Wert oder eine andere Teststatistik zurückgeben.
3. Verwenden Sie die replicateFunktion, um den Code ein paar Mal von oben auszuführen (ich beginne normalerweise mit ungefähr 100 Mal, um ein Gefühl dafür zu bekommen, wie lange es dauert und um den richtigen allgemeinen Bereich zu erhalten, dann bis zu 1.000 und manchmal 10.000 oder 100.000 für die Endwerte, die ich verwenden werde). Die Häufigkeit, mit der Sie die Nullhypothese zurückgewiesen haben, ist die Potenz.
4. Wiederholen Sie die obigen Schritte für weitere Bedingungen.

Hier ist ein einfaches Beispiel mit ordinaler Regression:

library(rms)

tmpfun <- function(n, beta0, beta1, beta2) {
    x <- runif(n, 0, 10)
    eta1 <- beta0 + beta1*x
    eta2 <- eta1 + beta2
    p1 <- exp(eta1)/(1+exp(eta1))
    p2 <- exp(eta2)/(1+exp(eta2))
    tmp <- runif(n)
    y <- (tmp < p1) + (tmp < p2)
    fit <- lrm(y~x)
    fit$stats[5]
}

out <- replicate(1000, tmpfun(100, -1/2, 1/4, 1/4))
mean( out < 0.05 )

Ich würde Snow's Antwort um eine weitere Sache ergänzen (und dies gilt für jede Leistungsanalyse mittels Simulation) - achten Sie darauf, ob Sie nach einem Test mit einem oder zwei Schwänzen suchen. Gängige Programme wie G * Power verwenden standardmäßig den 1-Tail-Test. Wenn Sie versuchen, festzustellen, ob Ihre Simulationen mit ihnen übereinstimmen (immer eine gute Idee, wenn Sie lernen, wie dies funktioniert), sollten Sie dies zuerst überprüfen.

Damit Snow einen 1-Tail-Test durchführt, würde ich den Funktionseingaben einen Parameter namens "tail" hinzufügen und der Funktion selbst etwas Ähnliches hinzufügen:

 #two-tail test
  if (tail==2) fit$stats[5]

  #one-tail test
  if (tail==1){
    if (fit$coefficients[5]>0) {
          fit$stats[5]/2
    } else 1

Die 1-tailed-Version prüft grundsätzlich, ob der Koeffizient positiv ist, und halbiert dann den p-Wert.

Neben dem hervorragenden Beispiel von Snow können Sie meines Erachtens auch eine Leistungssimulation durchführen, indem Sie aus einem vorhandenen Datensatz ein Resampling durchführen, das Ihren Effekt hat. Nicht ganz ein Bootstrap, da Sie nicht das gleiche n mit Ersatz abtasten , sondern die gleiche Idee.

Hier ist ein Beispiel: Ich führte ein kleines Selbstexperiment durch, das eine positive Punktschätzung ergab, bei der ordinalen logistischen Regression jedoch nicht annähernd statistisch signifikant war, da es sich um ein kleines Experiment handelte . Mit dieser Punkt-Schätzung, wie groß ein n brauche ich? Für verschiedene mögliche n habe ich oft einen Datensatz generiert und die ordinale logistische Regression durchgeführt und festgestellt, wie klein der p- Wert war:

library(boot)
library(rms)
npt <- read.csv("http://www.gwern.net/docs/nootropics/2013-gwern-noopept.csv")
newNoopeptPower <- function(dt, indices) {
    d <- dt[sample(nrow(dt), n, replace=TRUE), ] # new dataset, possibly larger than the original
    lmodel <- lrm(MP ~ Noopept + Magtein, data = d)
    return(anova(lmodel)[7])
}
alpha <- 0.05
for (n in seq(from = 300, to = 600, by = 30)) {
   bs <- boot(data=npt, statistic=newNoopeptPower, R=10000, parallel="multicore", ncpus=4)
   print(c(n, sum(bs$t<=alpha)/length(bs$t)))
}

Mit der Ausgabe (für mich):

[1] 300.0000   0.1823
[1] 330.0000   0.1925
[1] 360.0000   0.2083
[1] 390.0000   0.2143
[1] 420.0000   0.2318
[1] 450.0000   0.2462
[1] 480.000   0.258
[1] 510.0000   0.2825
[1] 540.0000   0.2855
[1] 570.0000   0.3184
[1] 600.0000   0.3175

In diesem Fall bei n betrug die Leistung = 600 32%. Nicht sehr ermutigend.

(Wenn mein Simulationsansatz falsch ist, sagen Sie es mir bitte. Ich werde ein paar medizinische Artikel über die Leistungssimulation zur Planung klinischer Studien veröffentlichen, bin mir jedoch nicht sicher, ob ich sie genau implementieren soll.)

In Bezug auf die erste Simulation (vorgeschlagen von Snow; /stats//a/22410/231675 ):

Ich bin mir immer noch nicht sicher, wie die Simulation mit mehr (genauer gesagt drei) unabhängigen Variablen aussehen soll. Ich verstehe, dass ich sie alle in das eta1-Teil aufnehmen sollte, z. B. eta1 <- beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3 '' (wie oben erwähnt). Ich weiß aber nicht, wie ich den Rest der Parameter in der Funktion einstellen soll. Könnte mir jemand dabei helfen?