R Power против кода силы сигнала

Question

Я надеюсь, что кто-то со свежими глазами сможет помочь мне здесь! Я пытаюсь раскрыть силу эксперимента, поэтому сделал следующее:

height <- seq(0, 151)
social <- rpois(length(height), 9 + 0.2 * (height))
m2 <- glm(score ~ height, family = poisson)
summary(m2)
m3 <- update(m2, ~. - height)
anova(m2, m3, test = "Chi")
test.results <- anova(m2, m3, test = "Chi")
names(test.results)
test.results$"P(>|Chi|)"
test.results$"P(>|Chi|)"[2]
get.p.value <- function(slope) {
    social <- rpois(length(height), 9 + slope * (height))
    m2 <- glm(score ~ height, family = poisson)
    m3 <- update(m2, ~. - r.hand)
    anova(m2, m3, test = "Chi")$"P(>|Chi|)"[2]
}
p.vals <- numeric(1000)
for (i in 1000) {
    p.vals[-0.5] <- get.p.value(-0.5)
}
p.vals
power.of.test <- length(p.vals[p.vals < 0.05])/length(p.vals)
power.of.test
slope.line <- seq(-0.2, -1.1, -0.1)
p.vals <- numeric(100)
power.of.test <- numeric(10)
for (j in 1:10) {
    for (i in 1:100) p.vals[i] <- get.p.value(slope.line[j])
    power.of.test[j] <- length(p.vals[p.vals < 0.05])/length(p.vals)
}
plot(slope.line, power.of.test)

Однако, это производит:

In rpois(length(height), 9 + slope * (height)) : NAs produced

Я, очевидно, где-то совершил глупую ошибку и провел весь день, перепечатывая ее, чтобы убедиться, что я не пропускаю скобки и т. Д., Но все, кажется, в порядке. У меня такое ощущение, что это как-то связано со значением 9 и уклоном, которое я получил из glm, но это может быть неправильно? Заранее спасибо.

Answer 1

Указанная вами ошибка находится в get.p.value(-0.5).Просмотр функции показывает, что вы присваиваете это (-0.5) значение rpois, которое не может принимать отрицательное число в качестве аргумента.

Я действительно не знаю, что вы хотите сделать, но вы можете запуститьваш код с get.p.value(0.5) вместо.

Удачи!

Answer 2

Несколько замечаний:

В функции glm у вас есть «оценка», но она не определена. Вы имеете в виду "социальный"? В цикле for у вас есть for(i in 1000), я думаю, вы имеете в виду for(i in 1:1000). Также

 p.vals[-0.5] <- get.p.value(-0.5)

не выполняет итерации по i, а p.vals[] обозначает индексирование ... что не имеет смысла, когда вы говорите о -5-м элементе p.vals.

Наконец строка

test.results$"P(>|Chi|)"[2]

дает мне результат 3,16 * 10 ^ -107, который на самом деле очень близок к 0. Это имеет смысл для меня, потому что я думаю, что вы спрашиваете код «удаляю ли я всю информацию из модели, какое значение иметь какую-либо информацию? "

Я мог бы быть здесь вне базы ... но, надеюсь, это поможет вам!

Answer 3

Попробуйте: сразу после for (j in 1:10) { введите browser() ... затем перед отправкой кода в командной строке наберите: options(warn=2)

Поэтому, когда вы получаете сообщение об ошибке или предупреждение (что, как я думаю, это сообщение), вы можете проверить любую из переменных, чтобы получить их текущие значения. Значение по умолчанию для warn равно 1, и вы должны сбросить его, когда закончите.

(Я подозреваю, что вместе с вами вы обнаружите, что 9 - (- some_slope) * some_height дает отрицательное значение для вызова rpois. (Ожидаемое значение для rpois действительно должно быть положительным , верно?)