0

I am trying to create 95% contour in R. I was advised to be more specific about my question in my recent post and hence posting it as a different question. I used an existing code from another post (R - How to find points within specific Contour) but am not sure why for levels I get an 'NA'. This is causing an error in the ls and point.in.polygon used later

library(MASS)
set.seed(1234)

p1 <- .136
p2 <- .069
nn <- 60

Y <- structure(list(Y1 = c(27L, 8L, 2L, 36L, 5L, 10L, 14L, 12L, 12L, 22L, 11L, 24L, 13L, 6L, 2L, 10L, 12L, 11L, 16L, 2L), Y2 = c(8L, 4L, 2L, 24L, 4L, 1L, 7L, 8L, 6L, 11L, 7L, 7L, 6L, 3L, 1L, 1L, 12L, 8L, 11L, 0L)), .Names = c("Y1", "Y2"), class = "data.frame", row.names = c(NA, -20L)) 

Y <- as.matrix(Y)
xx <- ifelse(Y==0,Y+.5,Y)
nnn <- ifelse(Y==0,nn+.5,nn)

xx <- as.matrix(xx)
Y1<-xx/nnn # estimates of p

sigma2 <- matrix(c(var(Y1[,1]),cov(Y1[,1],Y1[,2]),cov(Y1[,1],Y1[,2]),var(Y1[,2])),2,2)

rho <- sigma2[1,2]/sqrt(sigma2[1,1]*sigma2[2,2])

rate<-Y1[2,] # change for each site
print(rate)
#         Y1         Y2 
# 0.13333333 0.06666667 

rate1 <- rate/(1-rate)
rate2 <- log(rate1)

Sigma11 <- (1/(rate[1]*(1-rate[1]))^2)*sigma2[1,1]
Sigma22 <- (1/(rate[2]*(1-rate[2]))^2)*sigma2[2,2]
Sigma12 <- (1/((rate[1]*(1-rate[1]))*(rate[2]*(1-rate[2]))))*sigma2[1,2]

Sigma2 <- matrix(c(Sigma11,Sigma12,Sigma12,Sigma22),2,2)

print(Sigma2)
#          [,1]     [,2]
# [1,] 1.604805 1.503483
# [2,] 1.503483 2.102154

rate3 <- mvrnorm(1000, mu=c(rate2[1],rate2[2]), Sigma2)

x <- exp(rate3[,1])/(1+exp(rate3[,1]))
y <- exp(rate3[,2])/(1+exp(rate3[,2]))

## Points within polygons
library(MASS)
dens <- kde2d(x, y, n=100) 
image(dens)

prob <- c(0.975, 0.025)
dx <- diff(dens$x[1:2])
dy <- diff(dens$y[1:2])
sz <- sort(dens$z)

c1 <- cumsum(sz) * dx * dy 
levels <- sapply(prob, function(x) { 
  approx(c1, sz, xout = 1 - x)$y
})

print(levels)
# [1] 0.2480976        NA**   # <- Why NA?
#plot(p1,p2)
#contour(dens, levels=levels, labels=prob, add=T)
ls <- contourLines(dens, level=levels)

library(sp)
inner <- point.in.polygon(p1, p2, ls[[2]]$x, ls[[2]]$y) # whether points in inner ellipse
out <- point.in.polygon(p1, p2, ls[[1]]$x, ls[[1]]$y) # whether points in outter ellipse

within<-(inner+out)

print(within)
# [1] 0
Community
  • 1
  • 1
user1560215
  • 227
  • 1
  • 13
  • > dput(Y) structure(list(Y1 = c(27L, 8L, 2L, 36L, 5L, 10L, 14L, 12L, 12L, 22L, 11L, 24L, 13L, 6L, 2L, 10L, 12L, 11L, 16L, 2L), Y2 = c(8L, 4L, 2L, 24L, 4L, 1L, 7L, 8L, 6L, 11L, 7L, 7L, 6L, 3L, 1L, 1L, 12L, 8L, 11L, 0L)), .Names = c("Y1", "Y2"), class = "data.frame", row.names = c(NA, -20L)) – user1560215 Aug 09 '15 at 00:36
  • 1
    Create an [MCVE](http://stackoverflow.com/help/mcve) with emphasis on the M, and I bet you will have already solved your own problem. – A. Webb Aug 09 '15 at 01:54
  • It is the approx() function creating the problem which is meant for interpolation. When I feed 1-x with x be 0.025, the interpolation is outside of the range of ci, which will inevitably give us a NA as the result. For example, approx(1:100, 100:1, xout = 1000)$y will give NA. Other models as alternatives could be tried or an alternative extrapolation function – user1560215 Sep 10 '15 at 14:06

0 Answers0