[R] Need help with kde2d function in R

Thu Aug 6 02:25:21 CEST 2015

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]]