Hi All:
I am getting some warnings trying to run the following code and I am not
sure what they mean.
> set.seed(111)
> k<-63
> n<-60
> x<-NULL
> p<-rbeta(k,3,3)# so that the mean nausea rate is alpha/(alpha+beta)
> #min<-10
> #max<-60
> #n<-as.integer(runif(k,min,max))
> for(i in 1:k)
+ x<-cbind(x,rbinom(300,n,p[i]))
> x<-t(x)
> rate<-t(t(x)/n)
> se_rate<-sqrt(rate*(1-rate)/n)
>
>
>
> # Exact Confidence Interval
>
> l_cl_exact<-qbeta(.025,x,n-x+1)
There were 50 or more warnings (use warnings() to see the first 50)
> u_cl_exact<-qbeta(.975,x+1,n-x)
>
> for (i in 1:63){
+
+ for (j in 1:300)
+ {
+
+ if (x[i,j]==0)
+ {
+ l_cl_exact[i,j]<-0
+ u_cl_exact[i,j]<-u_cl_exact[i,j]
+ }
+ else if (x[i,j]==n[i])
+ {
+ l_cl_exact[i,j]<-l_cl_exact[i,j]
+ u_cl_exact[i,j]<-1
+ }
+ else
+ l_cl_exact[i,j]<-l_cl_exact[i,j]
+ u_cl_exact[i,j]<-u_cl_exact[i,j]
+
+ #print(c(i,j))
+
+ }
+ }
Error in if (x[i, j] == n[i]) { : missing value where TRUE/FALSE needed
>
>
>
>
>
>
> rate_t<-t(rate)
>
> l_cl_exact_t<-t(l_cl_exact)
> u_cl_exact_t<-t(u_cl_exact)
>
> coverage<-matrix(data=0,nrow=300,ncol=63)
> int_width<-matrix(data=0,nrow=300,ncol=63)
> sqr_err<-matrix(data=0,nrow=300,ncol=63)
> abs_bias<-matrix(data=0,nrow=300,ncol=63)
>
> for (i in 1:300){
+
+ for (j in 1:63)
+ {
+
+
+ if ((p[j]>l_cl_exact_t[i,j]) & (p[j]<u_cl_exact_t[i,j]))
+ {
+ coverage[i,j]<-1
+ }
+ else{
+ coverage[i,j]<-0
+
+ }
+
+ int_width[i,j]<-u_cl_exact_t[i,j]-l_cl_exact_t[i,j]
+
+ sqr_err[i,j]<-((rate_t[i,j]-p[j])^2)/0.5
+
+ abs_bias[i,j]<-abs(rate_t[i,j]-p[j])
+
+ #print(c(i,j))
+
+ }
+ }
> avg_cov<-apply(coverage,1,mean)
> cov_prob<-mean(avg_cov)
> 1-cov_prob# Non-Coverage Probability
[1] 0.03433862
>
> avg_int<-apply(int_width,1,mean)
> mean(avg_int)# Average Interval width
[1] 0.2408541
>
> avg_sqr_err<-apply(sqr_err,1,mean)
> mean(avg_sqr_err)# Average Squared Error Loss
[1] 0.007104948
>
> avg_abs_bias<-apply(abs_bias,1,mean)
> mean(abs_bias)
[1] 0.04738646
> warnings()
Warning messages:
1: In qbeta(0.025, x, n - x + 1) : pbeta_raw() -> bratio() gave error code 6
2: In qbeta(0.025, x, n - x + 1) : pbeta_raw() -> bratio() gave error code 6
3: In qbeta(0.025, x, n - x + 1) : pbeta_raw() -> bratio() gave error code 6
4: In qbeta(0.025, x, n - x + 1) : pbeta_raw() -> bratio() gave error code 6
5: In qbeta(0.025, x, n - x + 1) : pbeta_raw() -> bratio() gave error code 6
Thanks
Anamika
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
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.
