rchi-squaredbonferroni

How to calculate the Bonferroni Lower and Upper limits in R?


With the following data, I am trying to calculate the Chi Square and Bonferroni lower and upper Confidence intervals. The column "Data_No" identifies the dataset (as calculations needs to be done separately for each dataset).

Data_No    Area    Observed
   1        3353    31
   1        2297    2
   1        1590    15
   1        1087    16
   1        817     2
   1        847     10
   1        1014    28
   1        872     29
   1        1026    29
   1        1215    21
   2        3353    31
   2        2297    2
   2        1590    15
   3        1087    16
   3        817     2

The code I used is

        library(dplyr) 
        setwd("F:/GIS/July 2019/") 
        total_data <- read.csv("test.csv") 
        result_data <- NULL 
        for(i in unique(total_data$Data_No)){ 
        data <- total_data[which(total_data$Data_No == i),] data <- data %>%
        mutate(RelativeArea = Area/sum(Area), Expected = RelativeArea*sum(Observed), OminusE = Observed-Expected, O2 = OminusE^2, O2divE = O2/Expected, APU = Observed/sum(Observed), Alpha = 0.05/2*count(Data_No), 
Zvalue = qnorm(Alpha,lower.tail=FALSE), lower = APU-Zvalue*sqrt(APU*(1-APU)/sum(Observed)), upper = APU+Zvalue*sqrt(APU*(1-APU)/sum(Observed)))
result_data <- rbind(result_data,data) }
write.csv(result_data,file='final_result.csv')

And the error message I get is:

Error in UseMethod("summarise_") : no applicable method for 'summarise_' applied to an object of class "c('integer', 'numeric')"

The column that I am calling "Alpha" is the alpha value of 0.05/2k, where K is the number of categories - in my example, I have 10 categories ("Data_No" column) for the first dataset, so "Alpha" needs to be 0.05/20 = 0.0025, and it's corresponding Z value is 2.807. The second dataset has 3 categories (so 0.05/6) and the third has 2 categories (0.05/4) in my example table (Data_No" column). Using the values from the newly calculated "Alpha" column, I then need to calculate the ZValue column (Zvalue = qnorm(Alpha,lower.tail=FALSE)) which I then use to calculate the lower and upper confidence intervals.

From the above data, here are the results that I should get, but note that I have had to manually calculate Alpha column and Zvalue, rather than insert those calculations within the R code:

Data_No Area    Observed    RelativeArea    Alpha   Z value lower   upper
    1   3353    31          0.237           0.003   2.807   0.092   0.247
    1   2297    2           0.163           0.003   2.807   -0.011  0.033
    1   1590    15          0.113           0.003   2.807   0.025   0.139
    1   1087    16          0.077           0.003   2.807   0.029   0.146
    1   817     2           0.058           0.003   2.807   -0.011  0.033
    1   847     10          0.060           0.003   2.807   0.007   0.102
    1   1014    28          0.072           0.003   2.807   0.078   0.228
    1   872     29          0.062           0.003   2.807   0.083   0.234
    1   1026    29          0.073           0.003   2.807   0.083   0.234
    1   1215    21          0.086           0.003   2.807   0.049   0.181
    2   3353    31          0.463           0.008   2.394   0.481   0.811
    2   2297    2           0.317           0.008   2.394   -0.027  0.111
    2   1590    15          0.220           0.008   2.394   0.152   0.473
    3   1087    16          0.571           0.013   2.241   0.723   1.055
    3   817     2           0.429           0.013   2.241   -0.055  0.277

Please note that I only included some of the columns generated from the code.


Solution

  • # You need to check the closing bracket for lower c.f. sqrt value. Following code should work.
    
    data <- read.csv("test.csv") 
    data <- data %>% mutate(RelativeArea =
                              Area/sum(Area), Expected = RelativeArea*sum(Observed), OminusE =
                              Observed-Expected, O2 = OminusE^2, O2divE = O2/Expected, APU =
                              Observed/sum(Observed), lower =
                              APU-2.394*sqrt(APU*(1-APU)/sum(Observed)), upper =
                                               APU+2.394*sqrt(APU*(1-APU)/sum(Observed)))
    
    
    
    #Answer to follow-up question.
    #Sample Data
    Data_No   Area   Observed
    1         3353    31
    1         2297    2
    2         1590    15
    2         1087    16
    
    #Code to run
    total_data <- read.csv("test.csv")
    result_data <- NULL
    for(i in unique(total_data$Data_No)){
    data <- total_data[which(total_data$Data_No == i),]
    data <- data %>% mutate(RelativeArea =
                              Area/sum(Area), Expected = RelativeArea*sum(Observed), OminusE =
                              Observed-Expected, O2 = OminusE^2, O2divE = O2/Expected, APU =
                              Observed/sum(Observed), lower =
                              APU-2.394*sqrt(APU*(1-APU)/sum(Observed)), upper =
                                               APU+2.394*sqrt(APU*(1-APU)/sum(Observed)))
    
    result_data <- rbind(result_data,data)
    }
    
    write.csv(result_data,file='final_result.csv')