Problem with contour()

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Problem with contour()

Helmut Schütz
Dear all,

I can't get my head around how contour lines are drawn.
Working example(x and y are parameters of a certain test and z the
resulting power):

library(PowerTOST)
x    <- 0.90
y    <- 0.35
res  <- as.numeric(sampleN.TOST(theta0 = x, CV = y, design = "2x2x4",
                                 method = "central", details = FALSE,
                                 print = FALSE)[7:8])
mesh <- 28
ys   <- unique(sort(c(y, seq(y*.8, y*1.2, length.out = mesh))))
xs   <- unique(sort(c(x, seq(x*0.95, 1, length.out = mesh))))
z    <- matrix(nrow = length(ys), ncol = length(xs),
                dimnames = list(paste0("y.", signif(ys, 5)),
                                paste0("x.", signif(xs, 5))))
for (i in seq_along(ys)) {
   for (j in seq_along(xs)) {
     z[i, j] <- suppressMessages(
                  power.TOST(CV = ys[i], theta0 = xs[j], design = "2x2x4",
                             method = "central", n = res[1]))
   }
}
z    <- z[nrow(z):1, ncol(z):1]          # reverse rows & columns
z[paste0("y.", y), paste0("x.", x)] == res[2] # should be TRUE
contour(xs, ys, z, nlevels = 20, las = 1, labcex = 1,
         xlab = "x", ylab = "y", main = paste("n =", n))
abline(h = y, v = x, lty = 2)
points(x, y, col = "red", cex = 1.5)
text(x, y, labels = signif(z[paste0("y.", y), paste0("x.", x)], 6),
      col = "red", adj = c(-0.1, 1.5))

At x = 0.9 and y = 0.35 z = 0.8130092. Obviously this does not agree
with the contour-lines.
I'm sure that I screwed up - but where?

All the best,
Helmut

--
Ing. Helmut Schütz
BEBAC – Consultancy Services for
Bioequivalence and Bioavailability Studies
Neubaugasse 36/11
1070 Vienna, Austria
E [hidden email]
W https://bebac.at/
F https://forum.bebac.at/

______________________________________________
[hidden email] mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.