R Power vs Signal Strength code

I'm hoping someone with fresh eyes will be able to help me out here! I am attempting to discover the power of an experiment, so have done the following:

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 <-开发者_如何学C 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)

However, this produces:

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

I've obviously made a silly mistake somewhere and have spent all day retyping it to make sure I'm not missing a parenthesis, etc but everything seems to be in order. I have a feeling it's something to do with the 9 and slope value, which I obtained from a glm, but this could be wrong? Thanks in advance.

---

The error you pointed out is found in get.p.value(-0.5). Looking inside the function shows that your are giving this (-0.5) value to rpois, which cannot take negative number as argument.

I really do not know what you want to do, but you might run your code with get.p.value(0.5) instead.

Good luck!

---

Try this: right after for (j in 1:10) { put in browser() ... then before submitting the code at the command line type: options(warn=2)

So when you get an error or a warning (which I what I think that message is), you can inspect any of the variables to get their current values. The default value for warn is 1 and you should set it back when done.

(I suspect, along with you ,that you will find that 9-(-some_slope)*some_height gives a negative value for the rpois call. (The expected value to rpois does need to be positive, right?)

---

A few observations:

In the glm function you have "score" but that is not defined. Do you mean "social"? In the for loop you have for(i in 1000), I think you mean for(i in 1:1000). Also

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

doesn't iterate over i and p.vals[] indicates indexing...which doesn't make sense when you're talking about the -.5th element of p.vals.

Finally the line

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

gives me a result of 3.16 * 10^-107, which is real real close to 0. This makes sense to me, because I think you code is asking "if I remove all information from the model, what is the value of having any information?"

I could be off base here... but hopefully this helps you!