Hi,
I am trying to create a bivariate ellipse to see if the true values fall
within the 95% confidence ellipse. I am getting subscript out of bounds
error in R with following code. Not sure what is causing it. When I use the
kde2d function with n>=30 I get the error but not for n=25 and below.
library(MASS)
set.seed(1234)
#Set working directory
setwd("C://Tina/USB_Backup_042213/Paper II/MLN
Automation/csvs_equal_20s")
p1<- .136
p2<- .069
nn<-60
Y<-NULL
> Y <- read.csv(file=paste0("MVN",i,".csv"),
header=T)
>
> 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
>
>
> #print(Y1)
>
> sigma2<-
matrix(c(var(Y1[,1]),cov(Y1[,1],Y1[,2]),cov(Y1[,1],Y1[,2]),var(Y1[,2])),2,2)
> print(sigma2)
[,1] [,2]
[1,] 0.02142909 0.010810225
[2,] 0.01081023 0.008138709
>
> 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)
>
> #print(rate1)
>
> 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)
>
> rate3<-mvrnorm(1000, mu=c(rate2[1],rate2[2]), Sigma2)
>
> #print(rate3)
>
> 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=25)
> image(dens)
>
>
>
#filled.contour(dens,color.palette=colorRampPalette(c('white','blue','yellow','red','darkred')))
>
>
> 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
+ })
> #plot(p1,p2)
> #contour(dens, levels=levels, labels=prob, add=T)
> ls <- contourLines(dens, level=levels)
> print(ls)
[[1]]
[[1]]$level
[1] 0.2149866
[[1]]$x
[1] 0.004397130 0.004568786 0.020836478 0.032862816 0.040885850
0.049303040
[7] 0.077374571 0.095771788 0.113863291 0.139211514 0.127080458
0.139599553
[13] 0.150352011 0.186840732 0.201445193 0.223329452 0.239170012
0.245315258
[19] 0.259818172 0.296306893 0.332795613 0.369284333 0.405773053
0.415153288
[25] 0.442261774 0.472911758 0.478750494 0.515239214 0.536208449
0.551727935
[31] 0.588216655 0.624705375 0.646392119 0.624705375 0.618407028
0.607497171
[37] 0.624705375 0.661194096 0.697682816 0.734171536 0.750443060
0.770660257
[43] 0.780548310 0.807148977 0.823897778 0.815766169 0.807148977
0.770660257
[49] 0.745293482 0.750888227 0.734171536 0.733822194 0.734171536
0.770660257
[55] 0.780165097 0.807148977 0.821208006 0.807148977 0.770975720
0.770660257
[61] 0.734171536 0.712245087 0.697682816 0.693122349 0.693977032
0.697682816
[67] 0.704425811 0.719174523 0.700076998 0.706464923 0.711329752
0.697682816
[73] 0.661194096 0.624705375 0.588216655 0.577866122 0.588216655
0.589258511
[79] 0.588216655 0.551727935 0.515239214 0.478750494 0.469648285
0.442261774
[85] 0.405773053 0.395062033
[[1]]$y
[1] 0.1783278616 0.1786369552 0.2142104572 0.2497839592 0.2640206283
[6] 0.2853574612 0.3109011130 0.3209309632 0.3351406407 0.3565044651
[11] 0.3920779671 0.4276514691 0.4342808892 0.4573891723 0.4632249711
[16] 0.4823023723 0.4987984731 0.5343719750 0.5402153675 0.5529088973
[21] 0.5541833087 0.5511966122 0.5658454057 0.5699454770 0.5780502406
[26] 0.5699454770 0.5677285953 0.5582423238 0.5699454770 0.5766348870
[31] 0.5899083535 0.6007347908 0.6055189790 0.6249103884 0.6410924810
[36] 0.6766659830 0.6936330128 0.7051209093 0.7052210613 0.7104009871
[41] 0.7122394849 0.7175712767 0.7122394849 0.7002511608 0.6766659830
[46] 0.6410924810 0.6338424511 0.6130127829 0.6055189790 0.5699454770
[51] 0.5348306032 0.5343719750 0.5340918150 0.5053232008 0.4987984731
[56] 0.4728407406 0.4632249711 0.4535094747 0.4276514691 0.4274565392
[61] 0.4021096213 0.3920779671 0.3724826039 0.3565044651 0.3209309632
[66] 0.3068930025 0.2853574612 0.2497839592 0.2142104572 0.1786369552
[71] 0.1430634533 0.1356917738 0.1315167472 0.1273426210 0.1140927079
[76] 0.1074899513 0.0744396356 0.0719164493 0.0703319222 0.0509793010
[81] 0.0542138537 0.0417716905 0.0363429473 0.0207462171 0.0047801760
[86] 0.0007694453
[[2]]
[[2]]$level
[1] 0.2149866
[[2]]$x
[1] 0.2963069 0.2847163 0.2802502 0.2963069 0.3327956 0.3517197
0.3563875
[8] 0.3327956 0.2963069
[[2]]$y
[1] 0.5929265 0.6055190 0.6410925 0.6510478 0.6520754 0.6410925
0.6055190
[8] 0.5877331 0.5929265
[[3]]
[[3]]$level
[1] 0.2149866
[[3]]$x
[1] 0.8059980 0.8071490 0.8436377 0.8449198 0.8436377 0.8071490
0.7848487
[8] 0.7706603 0.7584362
[[3]]$y
[1] 0.8545335 0.8530500 0.8207634 0.8189600 0.8169512 0.8104161
0.8189600
[8] 0.8304354 0.8545335
> 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
Thanks
Anamika
[[alternative HTML version deleted]]