Chi-square Test

Loading...

Test of goodness of fit (Raw data)

One nominal variable

Note: Input values must be separated by tabs. Copy and paste from Excel/Numbers.

Your data needs to have the header (variable names) in the first row. Missing values should be indicated by a period (.) or NA.


                
                

Contingency table


                

Test result


                

Plot



R session info

              

Test of goodness of fit (Tabulated data)

One nominal variable

Note: Input values must be separated by tabs. Copy and paste from Excel/Numbers.

Your data needs to have the header (variable names) in the first row. Missing values should be indicated by a period (.) or NA.


                
                

Contingency table


                

Test result


                

Plot



R session info

              

Test of Independence (Raw data)

Two or more than two nominal variables

Note: Input values must be separated by tabs. Copy and paste from Excel/Numbers.

Your data needs to have the header (variable names) in the first row. Missing values should be indicated by a period (.) or NA.


                
                

Contingency table


                

Test result


                

Plot




R session info

              

Test of Independence (Tabulated data)

Two or more than two nominal variables

Note: Input values must be separated by tabs. Copy and paste from Excel/Numbers.

Your data needs to have the header (variable names) in the first row. Missing values should be indicated by a period (.) or NA.


                
                

Contingency table


                

Test result


                

Plot




R session info

              
Note

This web application is developed with Shiny.


List of Packages Used
library(shiny)
library(shinyAce)
library(pwr)
library(vcd)

Code

Source code for this application is based on "The handbook of Research in Foreign Language Learning and Teaching" (Takeuchi & Mizumoto, 2012).

The code for this web application is available at GitHub.

If you want to run this code on your computer (in a local R session), run the code below:
library(shiny)
runGitHub("chi","mizumot")

I referred to js-STAR for some parts of the codes. I would like to thank the authors of js-STAR, the very fast and excellent online software.


Citation in Publications

Mizumoto, A. (2015). Langtest (Version 1.0) [Web application]. Retrieved from http://langtest.jp


Article

Mizumoto, A., & Plonsky, L. (2015). R as a lingua franca: Advantages of using R for quantitative research in applied linguistics. Applied Linguistics, Advance online publication. doi:10.1093/applin/amv025


Recommended

To learn more about R, I suggest this excellent and free e-book (pdf), A Guide to Doing Statistics in Second Language Research Using R, written by Dr. Jenifer Larson-Hall.

Also, if you are a cool Mac user and want to use R with GUI, MacR is defenitely the way to go!


Author

Atsushi MIZUMOTO, Ph.D.
Associate Professor of Applied Linguistics
Faculty of Foreign Language Studies /
Graduate School of Foreign Language Education and Research,
Kansai University, Osaka, Japan



Code for "Chi-square Test"
by Atsushi Mizumoto

show with app
library(shiny)
library(shinyAce)
library(pwr)
library(vcd)



shinyServer(function(input, output) {



#----------------------------------------------------
# 1. Test of goodness of fit (Raw data)
#----------------------------------------------------

    data1 <- reactive({
        
        dat <- read.csv(text=input$text1, sep="", na.strings=c("","NA","."))
        
            x <- table(dat)
            n <- sum(x)
            x <- c(x, Sum=n)
            
            print(x)
        })
    
        output$data1.out <- renderPrint({
            data1()
        })
    
    
    
    
    
    test1 <- reactive({
        
        dat <- read.csv(text=input$text1, sep="", na.strings=c("","NA","."))
        
        
            x <- table(dat)
            chi <- chisq.test(x)
            print(chi)
        
            P0 <- rep(1/length(x), times = length(x))  # Expected ratio
            P1 <- x/sum(x)                             # Observed ratio
            w <- ES.w1(P0, P1)                         # Effect size w (Large=0.5, Medium=0.3, Small=0.1)
            
            cat("Effect size w =", round(w, 3), "\n")


            cat("\n", "---", "\n", "Multiple comparisons (p-value adjusted with Bonferroni method):", "\n", "\n")
            
            itemNum <- 0
            for (i in 1:length(chi$observed)) {
                for (j in 1:length(chi$observed)) {
                    if (i <= j) {next}
                    z <- (abs(chi$observed[[i]]-chi$observed[[j]])-1)/sqrt(chi$observed[[i]]+chi$observed[[j]])
                    p <- pnorm(z, lower.tail=FALSE)*2
                    p <- p.adjust(p, method = "bonferroni", n = length(x))
                    n <- chi$observed[[i]] + chi$observed[[j]]
                    p00 <- c(0.5, 0.5) # 母比率同等
                    p11 <- c(chi$observed[[i]]/n, chi$observed[[j]]/n) # 標本比率
                    p0  <- p00[1]
                    p1  <- p11[1]
                    ESg <- p1-p0 # 効果量g
                    cat(names(x)[j],"vs",names(x)[i],":", "z =", sprintf("%.3f",round(z,3)), ",", "p =", sprintf("%.3f",round(p,3)), ",", "Effect size g =", round(ESg,3), "\n")
                    itemNum <- itemNum + 1
                }
            }
    })
    
    output$test1.out <- renderPrint({
        test1()
    })
    
    
    
    
    
    makepPlot1 <- function(){
        
        dat <- read.csv(text=input$text1, sep="", na.strings=c("","NA","."))
        
            x <- table(dat)
            
            P0 <- rep(1/length(x), times = length(x))  # Expected ratio
            P1 <- x/sum(x)                             # Observed ratio
            
            par(mar=c(5,6,2,4))
            z <- matrix(c(P0, P1), nc=length(x), by=1)
            colnames(z) <- names(x)
            rownames(z) <- c("Expected", "Observed")
            barplot(t(z), hor=1, las=1, xlab="Percentage", col=gray.colors(length(x)))
            
            legend("bottomright",legend=colnames(z), fill=gray.colors(length(x)))
            
    }
    
    output$pPlot1 <- renderPlot({
        print(makepPlot1())
    })
    
    
    
    
    
    info1 <- reactive({
        info1 <- paste("This analysis was conducted with ", strsplit(R.version$version.string, " \\(")[[1]][1], ".", sep = "")
        info2 <- paste("It was executed on ", date(), ".", sep = "")
        cat(sprintf(info1), "\n")
        cat(sprintf(info2), "\n")
    })
    
    output$info1.out <- renderPrint({
        info1()
    })










#----------------------------------------------------
# 2. Test of goodness of fit (Tabulated data)
#----------------------------------------------------

    data2 <- reactive({
        
        dat <- read.csv(text=input$text2, sep="", na.strings=c("","NA","."))
        
            n <- sum(dat)
            x <- data.frame(dat, Sum = n)
            rownames(x) <- ""
            
            print(x)
        })
    
        output$data2.out <- renderPrint({
            data2()
        })
    
    
    
    
    
    test2 <- reactive({
        
        dat <- read.csv(text=input$text2, sep="", na.strings=c("","NA","."))
        
            x <- dat
            chi <- chisq.test(x)
            print(chi)
            
            P0 <- rep(1/length(x), times = length(x))  # Expected ratio
            P1 <- x/sum(x)                             # Observed ratio
            w <- ES.w1(P0, P1)                         # Effect size w (Large=0.5, Medium=0.3, Small=0.1)
        
        cat("Effect size w =", round(w, 3), "\n")
        
        
        cat("\n", "---", "\n", "Multiple comparisons (p-value adjusted with Bonferroni method):", "\n", "\n")
        
        itemNum <- 0
        for (i in 1:length(chi$observed)) {
            for (j in 1:length(chi$observed)) {
                if (i <= j) {next}
                z <- (abs(chi$observed[[i]]-chi$observed[[j]])-1)/sqrt(chi$observed[[i]]+chi$observed[[j]])
                p <- pnorm(z, lower.tail=FALSE)*2
                p <- p.adjust(p, method = "bonferroni", n = length(x))
                n <- chi$observed[[i]] + chi$observed[[j]]
                p00 <- c(0.5, 0.5) # 母比率同等
                p11 <- c(chi$observed[[i]]/n, chi$observed[[j]]/n) # 標本比率
                p0  <- p00[1]
                p1  <- p11[1]
                ESg <- p1-p0 # 効果量g
                cat(names(x)[j],"vs",names(x)[i],":", "z =", sprintf("%.3f",round(z,3)), ",", "p =", sprintf("%.3f",round(p,3)), ",", "Effect size g =", round(ESg,3), "\n")
                itemNum <- itemNum + 1
            }
        }
    })
    
    output$test2.out <- renderPrint({
        test2()
    })
    
    
    
    
    
    makepPlot2 <- function(){
        
        dat <- read.csv(text=input$text2, sep="", na.strings=c("","NA","."))
        
            x <- dat
        
            P0 <- rep(1/length(x), times = length(x))  # 期待比率
            P1 <- x/sum(x)               # 標本比率
        
            par(mar=c(5,6,2,4))
            z <- matrix(c(P0, P1), nc=length(x), by=1)
            colnames(z) <- names(x)
            rownames(z) <- c("Expected", "Observed")
            barplot(t(z), hor=1, las=1, xlab="Percentage", col=gray.colors(length(x)))
        
            legend("bottomright",legend=colnames(z), fill=gray.colors(length(x)))
    }
    
    output$pPlot2 <- renderPlot({
        print(makepPlot2())
    })
    
    
    
    
    
    info2 <- reactive({
        info1 <- paste("This analysis was conducted with ", strsplit(R.version$version.string, " \\(")[[1]][1], ".", sep = "")
        info2 <- paste("It was executed on ", date(), ".", sep = "")
        cat(sprintf(info1), "\n")
        cat(sprintf(info2), "\n")
    })
    
    output$info2.out <- renderPrint({
        info2()
    })










#----------------------------------------------------
# 3. Test of Independence (Raw data)
#----------------------------------------------------

    data3 <- reactive({
        
        dat <- read.csv(text=input$text3, sep="", na.strings=c("","NA","."))
        
            x <- table(dat)
            x <- addmargins(x)
        
            print(x)
        })
    
        output$data3.out <- renderPrint({
            data3()
        })
    
    
    
    
    
    test3 <- reactive({
        
        dat <- read.csv(text=input$text3, sep="", na.strings=c("","NA","."))
        
            x <- table(dat)
            
            a <- chisq.test(x, correct=F)           # Pearson's Chi-squared
            b <- chisq.test(x)                      # Yates
            c <- assocstats(x)                      # Likelihood Ratio
            
            
            aa <- data.frame(a[4], a[1], a[2], a[3])
            aa[1] <- c("Pearson's Chi-squared")
            row.names(aa) <- NULL
            
            bb <- data.frame(b[4], b[1], b[2], b[3])
            bb[1] <- c("Yates' Continuity Correction")
            row.names(bb) <- NULL
            
            cc <- data.frame(c[2])[1,]
            cc <- data.frame(c("Log-likelihood Ratio (G)"), cc[,1], cc[,2], cc[,3])
            names(cc) <- c("method", "statistic", "parameter", "p.value")
            
            
            # Fisher's Exact Test
            e <- try(fisher.test(x, workspace=3000000), silent = FALSE)   # try, fisher.test
            if (class(e) == "try-error") {
                dd <- data.frame(c(""), c(""), c(""), c("data too large"))  # if error
                dd[1] <- c("Fisher's Exact Test")
                names(dd) <- c("method", "statistic", "parameter", "p.value")
                row.names(dd) <- NULL
            } else {
                d <- fisher.test(x, workspace=3000000)    # no error
                if (nrow(x) * ncol(x) == 4) { # Only for the 2×2 table
                    dd <- data.frame(d[6], c(""), c(""), d[1])
                } else {
                    dd <- data.frame(d[3], c(""), c(""), d[1])
                }
                dd[1] <- c("Fisher's Exact Test")
                names(dd) <- c("method", "statistic", "parameter", "p.value")
                row.names(dd) <- NULL
            }
            
            
            res <- rbind(aa, bb, cc, dd)
            res[2] <- round(as.numeric(res[,2]), 4)
            res[2][4,] <- ""
            names(res) <- c("Test", "X-squared", "df", "p-value")
            print(res)
            
            
            cat("\n")
            cat("\n", "------------------------------------------------------", "\n",
            "Effect size:", "\n")
            
            # Cramer's V [95%CI]
            V <- sqrt(chisq.test(x, correct=F)$statistic[[1]]/(sum(x)*(min(nrow(x),ncol(x)-1))))
            fisherZ <- 0.5 * log((1 + V)/(1 - V))
            fisherZ.se <- 1/sqrt(sum(x) - 3) * qnorm(1 - ((1 - 0.95)/2))
            fisherZ.ci <- fisherZ + c(-fisherZ.se, fisherZ.se)
            ci.V <- (exp(2 * fisherZ.ci) - 1)/(1 + exp(2 * fisherZ.ci))
            
            cat("\n", "Cramer's V [95%CI] =", V, "[", ci.V[1], ",",ci.V[2],"]", "\n")

            # Odds Ratio
            if (nrow(x) * ncol(x) == 4) { # Print odds ratio only for the 2×2 table
                
                oddsrt <- (x[1,1]/x[1,2]) / (x[2,1]/x[2,2])
                oddsrt.log <- log((x[1,1]/x[1,2]) / (x[2,1]/x[2,2]))
                oddsrt.log.var <- 1/x[1,1] + 1/x[1,2] + 1/x[2,1] + 1/x[2,2]
                ci.lwrupr <- exp(oddsrt.log + qnorm(c(0.025,0.975)) * sqrt(oddsrt.log.var))
                cat("\n", "Odds Ratio [95%CI] =", oddsrt, "[", ci.lwrupr[1], ",",ci.lwrupr[2],"]", "\n")
                
            } else {
                NULL
            }
            
            
            
            cat("\n", "\n", "------------------------------------------------------", "\n", "Residual analysis:", "\n")
            res <- chisq.test(x)  # 検定結果を代入
            
            cat("\n", "[Expected values]", "\n")
            print(res$expected) # 期待値
            
            cat("\n", "[Standardized residuals]", "\n")
            print(res$residuals) # 標準化残差
            
            cat("\n", "[Adjusted standardized residuals]", "\n")
            print(res$residuals/sqrt(outer(1-rowSums(x)/sum(x), 1-colSums(x)/sum(x)))) # 調整済み標準化残差
            
            cat("\n", "[p-values of adjusted standardized residuals (two-tailed)]", "\n")
            print(round(2*(1-pnorm(abs(res$residuals/sqrt(outer(1-rowSums(x)/sum(x), 1-colSums(x)/sum(x)))))),3)) # 残差の調整後有意確率(両側確率)
            
            
            
            # 多重比較
            if (nrow(x) < 3) { # Not showing if the number of row is 2
                
                cat("\n")
                
            } else {
                
                cat("\n", "\n", "------------------------------------------------------", "\n", "Multiple comparisons (p-value adjusted with Bonferroni method):", "\n", "\n")
                
                levI <- nrow(x)
                levJ <- ncol(x)
                N <- sum(x)
                dosu <- as.vector(t(x))
                M <- min(c(levI, levJ))
                
                V.95CI <- function(x) {
                    V <- sqrt(chisq.test(x, correct=F)$statistic[[1]]/(sum(x)*(min(nrow(x),ncol(x)-1))))
                    fisherZ <- 0.5 * log((1 + V)/(1 - V))
                    fisherZ.se <- 1/sqrt(sum(x) - 3) * qnorm(1 - ((1 - 0.95)/2))
                    fisherZ.ci <- fisherZ + c(-fisherZ.se, fisherZ.se)
                    ci.V <- (exp(2 * fisherZ.ci) - 1)/(1 + exp(2 * fisherZ.ci))
                    return(ci.V)
                }
                
                # エラーメッセージが出る場合の例外処理
                e <- try(fisher.test(x, workspace=3000000), silent = FALSE)   # try でfisher.test関数を実行
                if (class(e) == "try-error") { # if error, without fisher.p
                    
                    a <- c()
                    stat <- c()
                    jiyu <- c()
                    pchi <- c()
                    v <- c()
                    v.lower <- c()
                    v.upper <- c()
                    for(i in 1:(levI-1))
                    {
                        for(k in (i+1):levI)
                        {
                            ds <- c()
                            for(j in 1:levJ)
                            {
                                ds <- c(ds, dosu[(i-1)*levJ+j])
                            }
                            for(j in 1:levJ)
                            {
                                ds <- c(ds, dosu[(k-1)*levJ+j])
                            }
                            
                            kaik <- c()
                            kaik <- chisq.test(matrix(ds, nr=2, by=1), correct=F)
                            stat <- c(stat,kaik$stat)
                            jiyu <- c(jiyu,kaik$para)
                            pchi <- c(pchi,kaik$p.va)
                            v <- c(v, sqrt(kaik$stat/(N*(M-1))))
                            ci.values <- V.95CI(matrix(ds, nr=2, by=1))
                            lower <- ci.values[1]
                            upper <- ci.values[2]
                            v.lower <- c(v.lower, lower)
                            v.upper <- c(v.upper, upper)
                        } }
                    padj <- c()
                    padj <- p.adjust(pchi, "bonferroni")
                    
                    kochi<-c(); aite<-c()
                    for(i in 1:(levI-1))
                    {
                        for(j in (i+1):levI)
                        {
                            kochi <- c(kochi, rownames(x)[i])
                            aite  <- c(aite,  rownames(x)[j])
                        }
                    }
                    
                    a <- data.frame(Comparisons=paste(kochi, aite, sep=" vs "), "X^2"=round(stat, 4),
                    df=round(jiyu, 4),
                    p=round(padj, 4),
                    Cramers.V=round(v, 4),
                    lwr.V=round(v.lower, 4),
                    upr.V=round(v.upper, 4))
                    rownames(a) <- cat("\n")
                    
                    print(a)
                    
                } else { # with fisher.p
                    
                    
                    a <- c()
                    stat <- c()
                    jiyu <- c()
                    pchi <- c()
                    fishp <- c()
                    v <- c()
                    v.lower <- c()
                    v.upper <- c()
                    for(i in 1:(levI-1))
                    {
                        for(k in (i+1):levI)
                        {
                            ds <- c()
                            for(j in 1:levJ)
                            {
                                ds <- c(ds, dosu[(i-1)*levJ+j])
                            }
                            for(j in 1:levJ)
                            {
                                ds <- c(ds, dosu[(k-1)*levJ+j])
                            }
                            
                            kaik <- c()
                            kaik <- chisq.test(matrix(ds, nr=2, by=1), correct=F)
                            stat <- c(stat,kaik$stat)
                            jiyu <- c(jiyu,kaik$para)
                            pchi <- c(pchi,kaik$p.va)
                            fishp <- c(fishp,fisher.test(matrix(ds, nr=2, by=1))$p.va)
                            v <- c(v, sqrt(kaik$stat/(N*(M-1))))
                            ci.values <- V.95CI(matrix(ds, nr=2, by=1))
                            lower <- ci.values[1]
                            upper <- ci.values[2]
                            v.lower <- c(v.lower, lower)
                            v.upper <- c(v.upper, upper)
                        } }
                    padj <- c()
                    padj <- p.adjust(pchi, "bonferroni")
                    fishpadj <- p.adjust(fishp, "bonferroni")
                    
                    kochi<-c(); aite<-c()
                    for(i in 1:(levI-1))
                    {
                        for(j in (i+1):levI)
                        {
                            kochi <- c(kochi, rownames(x)[i])
                            aite  <- c(aite,  rownames(x)[j])
                        }
                    }
                    
                    a <- data.frame(Comparisons=paste(kochi, aite, sep=" vs "), "X^2"=round(stat, 4),
                    df=round(jiyu, 4),
                    p=round(padj, 4),
                    Fisher.p=round(fishpadj, 4),
                    Cramers.V=round(v, 4),
                    lwr.V=round(v.lower, 4),
                    upr.V=round(v.upper, 4))
                    rownames(a) <- cat("\n")
                    
                    print(a)
            }
        }
    })
    
    output$test3.out <- renderPrint({
        test3()
    })
    
    
    
    
    
    makepPlot3 <- function(){
        
        dat <- read.csv(text=input$text3, sep="", na.strings=c("","NA","."))
        
            x <- table(dat)
            
            levI <- nrow(x) # 行の水準数
            levJ <- ncol(x) # 列の水準数
            dosu <- as.vector(t(x))
            
            # 標本比率の計算
            gokei <- c()
            bunbo <- c()
            for(i in 1:levI) # 各群の度数を集計
            {
                ds <- c()
                for(j in 1:levJ)
                {
                    ds <- c(ds, dosu[(i-1)*levJ+j])
                }
                gokei <- c(gokei, sum(ds))
                bunbo <- c(bunbo, rep(sum(ds), levJ))
            }
            hyohir <- dosu/bunbo # 群別の各値の比率
            
            zuhir <- c()
            for(i in levI:1) # 群i→群1と逆順に並べ替える
            {
                for(j in 1:levJ)
                {
                    zuhir <- c(zuhir, hyohir[(i-1)*levJ+j] )
                }
            }
            
            zubar <- matrix(c(zuhir), nc=levJ, by=1)
            rownames(zubar) <- rev(rownames(x))
            colnames(zubar) <- colnames(x)
            #zubar <- zubar[nrow(zubar):1,]
            
            # プロット
            par(mar=c(5,6,2,4))
            barplot(t(zubar), hor=1, las=1, xlab="Percentage", col=gray.colors(ncol(x)))
            legend("bottomright", legend=colnames(zubar), fill=gray.colors(ncol(x)))
    }
    
    output$pPlot3 <- renderPlot({
        print(makepPlot3())
    })
    
    
    
    
    
   makemPlot3 <- function(){
       
       dat <- read.csv(text=input$text3, sep="", na.strings=c("","NA","."))
       
           x <- table(dat)
           mosaic(x, gp = shading_max, main="Mosaic plot")
       
    }
    
    output$mPlot3 <- renderPlot({
        print(makemPlot3())
    })
    
    



    info3 <- reactive({
        info1 <- paste("This analysis was conducted with ", strsplit(R.version$version.string, " \\(")[[1]][1], ".", sep = "")
        info2 <- paste("It was executed on ", date(), ".", sep = "")
        cat(sprintf(info1), "\n")
        cat(sprintf(info2), "\n")
    })
    
    output$info3.out <- renderPrint({
        info3()
    })










#----------------------------------------------------
# 4. Test of Independence (Tabulated data)
#----------------------------------------------------

    data4 <- reactive({
        
        dat <- read.csv(text=input$text4, sep="", na.strings=c("","NA","."))
        
            x <- as.matrix(dat)
            x <- addmargins(x)
        
            print(x)
        })
    
        output$data4.out <- renderPrint({
            data4()
        })
    
    
    
    
    
    test4 <- reactive({
        
        dat <- read.csv(text=input$text4, sep="", na.strings=c("","NA","."))
        
            x <- as.matrix(dat)
            
            
            a <- chisq.test(x, correct=F)           # Pearson's Chi-squared
            b <- chisq.test(x)                      # Yates
            c <- assocstats(x)                      # Likelihood Ratio
            
            aa <- data.frame(a[4], a[1], a[2], a[3])
            aa[1] <- c("Pearson's Chi-squared")
            row.names(aa) <- NULL
            
            bb <- data.frame(b[4], b[1], b[2], b[3])
            bb[1] <- c("Yates' Continuity Correction")
            row.names(bb) <- NULL
            
            cc <- data.frame(c[2])[1,]
            cc <- data.frame(c("Log-likelihood Ratio (G)"), cc[,1], cc[,2], cc[,3])
            names(cc) <- c("method", "statistic", "parameter", "p.value")
            
            
            # Fisher's Exact Test
            e <- try(fisher.test(x, workspace=3000000), silent = FALSE)   # try, fisher.test
            if (class(e) == "try-error") {
                dd <- data.frame(c(""), c(""), c(""), c("data too large"))  # if error
                dd[1] <- c("Fisher's Exact Test")
                names(dd) <- c("method", "statistic", "parameter", "p.value")
                row.names(dd) <- NULL
            } else {
                d <- fisher.test(x, workspace=3000000)    # no error
                if (nrow(x) * ncol(x) == 4) { # Only for the 2×2 table
                    dd <- data.frame(d[6], c(""), c(""), d[1])
                } else {
                    dd <- data.frame(d[3], c(""), c(""), d[1])
                }
                dd[1] <- c("Fisher's Exact Test")
                names(dd) <- c("method", "statistic", "parameter", "p.value")
                row.names(dd) <- NULL
            }
            
            
            res <- rbind(aa, bb, cc, dd)
            res[2] <- round(as.numeric(res[,2]), 4)
            res[2][4,] <- ""
            names(res) <- c("Test", "X-squared", "df", "p-value")
            cat("\n")
            print(res)
            
            
            
            
            cat("\n")
            cat("\n", "------------------------------------------------------", "\n",
            "Effect size:", "\n")

            # Cramer's V [95%CI]
            V <- sqrt(chisq.test(x, correct=F)$statistic[[1]]/(sum(x)*(min(nrow(x),ncol(x)-1))))
            fisherZ <- 0.5 * log((1 + V)/(1 - V))
            fisherZ.se <- 1/sqrt(sum(x) - 3) * qnorm(1 - ((1 - 0.95)/2))
            fisherZ.ci <- fisherZ + c(-fisherZ.se, fisherZ.se)
            ci.V <- (exp(2 * fisherZ.ci) - 1)/(1 + exp(2 * fisherZ.ci))
            
            cat("\n", "Cramer's V [95%CI] =", V, "[", ci.V[1], ",",ci.V[2],"]", "\n")
            
            
            # Odds Ratio
            if (nrow(x) * ncol(x) == 4) { # Print odds ratio only for the 2×2 table
               
               oddsrt <- (x[1,1]/x[1,2]) / (x[2,1]/x[2,2])
               oddsrt.log <- log((x[1,1]/x[1,2]) / (x[2,1]/x[2,2]))
               oddsrt.log.var <- 1/x[1,1] + 1/x[1,2] + 1/x[2,1] + 1/x[2,2]
               ci.lwrupr <- exp(oddsrt.log + qnorm(c(0.025,0.975)) * sqrt(oddsrt.log.var))
               cat("\n", "Odds Ratio [95%CI] =", oddsrt, "[", ci.lwrupr[1], ",",ci.lwrupr[2],"]", "\n")
               
            } else {
                NULL
            }
            
            
            
            cat("\n", "\n", "------------------------------------------------------", "\n", "Residual analysis:", "\n")
            res <- chisq.test(x)  # 検定結果を代入
            
            cat("\n", "[Expected values]", "\n")
            print(res$expected) # 期待値
            
            cat("\n", "[Standardized residuals]", "\n")
            print(res$residuals) # 標準化残差
            
            cat("\n", "[Adjusted standardized residuals]", "\n")
            print(res$residuals/sqrt(outer(1-rowSums(x)/sum(x), 1-colSums(x)/sum(x)))) # 調整済み標準化残差
            
            cat("\n", "[p-values of adjusted standardized residuals (two-tailed)]", "\n")
            print(round(2*(1-pnorm(abs(res$residuals/sqrt(outer(1-rowSums(x)/sum(x), 1-colSums(x)/sum(x)))))),3)) # 残差の調整後有意確率(両側確率)
            
            
            
            # 多重比較
            if (nrow(x) < 3) { # Not showing if the number of row is 2
                
                cat("\n")
                
            } else {
                
            cat("\n", "\n", "------------------------------------------------------", "\n", "Multiple comparisons (p-value adjusted with Bonferroni method):", "\n", "\n")
            
            levI <- nrow(x)
            levJ <- ncol(x)
            N <- sum(x)
            dosu <- as.vector(t(x))
            M <- min(c(levI, levJ))
            
            V.95CI <- function(x) {
                V <- sqrt(chisq.test(x, correct=F)$statistic[[1]]/(sum(x)*(min(nrow(x),ncol(x)-1))))
                fisherZ <- 0.5 * log((1 + V)/(1 - V))
                fisherZ.se <- 1/sqrt(sum(x) - 3) * qnorm(1 - ((1 - 0.95)/2))
                fisherZ.ci <- fisherZ + c(-fisherZ.se, fisherZ.se)
                ci.V <- (exp(2 * fisherZ.ci) - 1)/(1 + exp(2 * fisherZ.ci))
                return(ci.V)
            }
            
            # エラーメッセージが出る場合の例外処理
            e <- try(fisher.test(x, workspace=3000000), silent = FALSE)   # try でfisher.test関数を実行
            if (class(e) == "try-error") { # if error, without fisher.p
                
                a <- c()
                stat <- c()
                jiyu <- c()
                pchi <- c()
                v <- c()
                v.lower <- c()
                v.upper <- c()
                for(i in 1:(levI-1))
                {
                    for(k in (i+1):levI)
                    {
                        ds <- c()
                        for(j in 1:levJ)
                        {
                            ds <- c(ds, dosu[(i-1)*levJ+j])
                        }
                        for(j in 1:levJ)
                        {
                            ds <- c(ds, dosu[(k-1)*levJ+j])
                        }
                        
                        kaik <- c()
                        kaik <- chisq.test(matrix(ds, nr=2, by=1), correct=F)
                        stat <- c(stat,kaik$stat)
                        jiyu <- c(jiyu,kaik$para)
                        pchi <- c(pchi,kaik$p.va)
                        v <- c(v, sqrt(kaik$stat/(N*(M-1))))
                        ci.values <- V.95CI(matrix(ds, nr=2, by=1))
                        lower <- ci.values[1]
                        upper <- ci.values[2]
                        v.lower <- c(v.lower, lower)
                        v.upper <- c(v.upper, upper)
                    } }
                padj <- c()
                padj <- p.adjust(pchi, "bonferroni")
                
                kochi<-c(); aite<-c()
                for(i in 1:(levI-1))
                {
                    for(j in (i+1):levI)
                    {
                        kochi <- c(kochi, rownames(x)[i])
                        aite  <- c(aite,  rownames(x)[j])
                    }
                }
                
                a <- data.frame(Comparisons=paste(kochi, aite, sep=" vs "), "X^2"=round(stat, 4),
                df=round(jiyu, 4),
                p=round(padj, 4),
                Cramers.V=round(v, 4),
                lwr.V=round(v.lower, 4),
                upr.V=round(v.upper, 4))
                rownames(a) <- cat("\n")
                
                print(a)
                
            } else { # with fisher.p
                
                
                a <- c()
                stat <- c()
                jiyu <- c()
                pchi <- c()
                fishp <- c()
                v <- c()
                v.lower <- c()
                v.upper <- c()
                for(i in 1:(levI-1))
                {
                    for(k in (i+1):levI)
                    {
                        ds <- c()
                        for(j in 1:levJ)
                        {
                            ds <- c(ds, dosu[(i-1)*levJ+j])
                        }
                        for(j in 1:levJ)
                        {
                            ds <- c(ds, dosu[(k-1)*levJ+j])
                        }
                        
                        kaik <- c()
                        kaik <- chisq.test(matrix(ds, nr=2, by=1), correct=F)
                        stat <- c(stat,kaik$stat)
                        jiyu <- c(jiyu,kaik$para)
                        pchi <- c(pchi,kaik$p.va)
                        fishp <- c(fishp,fisher.test(matrix(ds, nr=2, by=1))$p.va)
                        v <- c(v, sqrt(kaik$stat/(N*(M-1))))
                        ci.values <- V.95CI(matrix(ds, nr=2, by=1))
                        lower <- ci.values[1]
                        upper <- ci.values[2]
                        v.lower <- c(v.lower, lower)
                        v.upper <- c(v.upper, upper)
                    } }
                padj <- c()
                padj <- p.adjust(pchi, "bonferroni")
                fishpadj <- p.adjust(fishp, "bonferroni")
                
                kochi<-c(); aite<-c()
                for(i in 1:(levI-1))
                {
                    for(j in (i+1):levI)
                    {
                        kochi <- c(kochi, rownames(x)[i])
                        aite  <- c(aite,  rownames(x)[j])
                    }
                }
                
                a <- data.frame(Comparisons=paste(kochi, aite, sep=" vs "), "X^2"=round(stat, 4),
                df=round(jiyu, 4),
                p=round(padj, 4),
                Fisher.p=round(fishpadj, 4),
                Cramers.V=round(v, 4),
                lwr.V=round(v.lower, 4),
                upr.V=round(v.upper, 4))
                rownames(a) <- cat("\n")
                
                print(a)
            }
            }
    })
    
    output$test4.out <- renderPrint({
        test4()
    })
    
    
    
    
    
    makepPlot4 <- function(){
        
        dat <- read.csv(text=input$text4, sep="", na.strings=c("","NA","."))
        
           x <- as.matrix(dat)
            
            levI <- nrow(x) # 行の水準数
            levJ <- ncol(x) # 列の水準数
            dosu <- as.vector(t(x))
            
            # 標本比率の計算
            gokei <- c()
            bunbo <- c()
            for(i in 1:levI) # 各群の度数を集計
            {
                ds <- c()
                for(j in 1:levJ)
                {
                    ds <- c(ds, dosu[(i-1)*levJ+j])
                }
                gokei <- c(gokei, sum(ds))
                bunbo <- c(bunbo, rep(sum(ds), levJ))
            }
            hyohir <- dosu/bunbo # 群別の各値の比率
            
            zuhir <- c()
            for(i in levI:1) # 群i→群1と逆順に並べ替える
            {
                for(j in 1:levJ)
                {
                    zuhir <- c(zuhir, hyohir[(i-1)*levJ+j] )
                }
            }
            
            zubar <- matrix(c(zuhir), nc=levJ, by=1)
            rownames(zubar) <- rev(rownames(x))
            colnames(zubar) <- colnames(x)
            #zubar <- zubar[nrow(zubar):1,]
            
            # プロット
            par(mar=c(5,6,2,4))
            barplot(t(zubar), hor=1, las=1, xlab="Percentage", col=gray.colors(ncol(x)))
            legend("bottomright", legend=colnames(zubar), fill=gray.colors(ncol(x)))
    
    }
    
    output$pPlot4 <- renderPlot({
        print(makepPlot4())
    })
    
    
    
    
    
    makemPlot4 <- function(){
       
       dat <- read.csv(text=input$text4, sep="", na.strings=c("","NA","."))
       
            x <- as.matrix(dat)
            mosaic(x, gp = shading_max, main="Mosaic plot")
       
    }
    
    output$mPlot4 <- renderPlot({
        print(makemPlot4())
    })
    
    



    info4 <- reactive({
        info1 <- paste("This analysis was conducted with ", strsplit(R.version$version.string, " \\(")[[1]][1], ".", sep = "")
        info2 <- paste("It was executed on ", date(), ".", sep = "")
        cat(sprintf(info1), "\n")
        cat(sprintf(info2), "\n")
    })
    
    output$info4.out <- renderPrint({
        info4()
    })






})
library(shiny)
library(shinyAce)


shinyUI(bootstrapPage(


    headerPanel("Chi-square Test"),


########## Adding loading message #########

tags$head(tags$style(type="text/css", "
#loadmessage {
position: fixed;
top: 0px;
left: 0px;
width: 100%;
padding: 10px 0px 10px 0px;
text-align: center;
font-weight: bold;
font-size: 100%;
color: #000000;
background-color: #CCFF66;
z-index: 105;
}
")),

conditionalPanel(condition="$('html').hasClass('shiny-busy')",
tags$div("Loading...",id="loadmessage")),

########## Added up untill here ##########



    mainPanel(
        tabsetPanel(position = "left", selected = "Test of Independence (Tabulated data)",

        tabPanel("Test of goodness of fit (Raw data)",

            h2("Test of goodness of fit (Raw data)"),

            h4("One nominal variable"),

            p('Note: Input values must be separated by tabs. Copy and paste from Excel/Numbers.'),

            p(HTML("<b><div style='background-color:#FADDF2;border:1px solid black;'>Your data needs to have the header (variable names) in the first row. Missing values should be indicated by a period (.) or NA.</div></b>")),

            aceEditor("text1", value="L1\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nJapanese\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nThai\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese\nChinese", mode="r", theme="cobalt"),

            br(),

            h3("Contingency table"),
            verbatimTextOutput("data1.out"),

            br(),

            h3("Test result"),
            verbatimTextOutput("test1.out"),

            br(),

            h3("Plot"),

            plotOutput("pPlot1"),

            br(),
            br(),

            strong('R session info'),
            verbatimTextOutput("info1.out")
            ),









        tabPanel("Test of goodness of fit (Tabulated data)",

            h2("Test of goodness of fit (Tabulated data)"),

            h4("One nominal variable"),

            p('Note: Input values must be separated by tabs. Copy and paste from Excel/Numbers.'),

            p(HTML("<b><div style='background-color:#FADDF2;border:1px solid black;'>Your data needs to have the header (variable names) in the first row. Missing values should be indicated by a period (.) or NA.</div></b>")),

            aceEditor("text2", value="Japanese\tThai\tChinese\n18\t24\t48", mode="r", theme="cobalt"),

            br(),

            h3("Contingency table"),
            verbatimTextOutput("data2.out"),

            br(),

            h3("Test result"),
            verbatimTextOutput("test2.out"),

            br(),

            h3("Plot"),

            plotOutput("pPlot2"),

            br(),
            br(),

            strong('R session info'),
            verbatimTextOutput("info2.out")
            ),










        tabPanel("Test of Independence (Raw data)",

            h2("Test of Independence (Raw data)"),

            h4("Two or more than two nominal variables"),

            p('Note: Input values must be separated by tabs. Copy and paste from Excel/Numbers.'),

            p(HTML("<b><div style='background-color:#FADDF2;border:1px solid black;'>Your data needs to have the header (variable names) in the first row. Missing values should be indicated by a period (.) or NA.</div></b>")),

            aceEditor("text3", value="Sex\tEffect\nM\tNo\nW\tNo\nW\tNo\nM\tNo\nM\tYes\nM\tYes\nM\tYes\nM\tNo\nW\tYes\nM\tNo\nW\tYes\nM\tNo\nM\tYes\nM\tNo\nM\tNo\nM\tYes\nW\tYes\nW\tYes\nW\tYes\nW\tYes\nW\tYes\nM\tYes\nM\tNo\nM\tNo\nM\tYes\nM\tYes\nW\tYes\nM\tNo\nM\tYes\nW\tYes\nM\tNo\nM\tNo\nW\tYes\nW\tYes\nW\tYes\nW\tYes\nM\tNo\nW\tNo\nW\tYes\nM\tYes\nW\tYes\nM\tNo\nM\tYes\nW\tYes\nM\tYes\nW\tYes\nM\tYes\nM\tNo\nM\tNo\nW\tNo\nW\tNo\nM\tYes\nW\tNo\nM\tYes\nW\tYes\nW\tYes\nM\tNo\nM\tNo\nM\tYes\nW\tYes\nM\tNo\nW\tYes\nW\tYes\nM\tYes\nW\tNo\nW\tYes\nM\tNo\nW\tYes\nW\tNo\nM\tYes",mode="r", theme="cobalt"),

            br(),

            h3("Contingency table"),
            verbatimTextOutput("data3.out"),

            br(),

            h3("Test result"),
            verbatimTextOutput("test3.out"),

            br(),

            h3("Plot"),

            plotOutput("pPlot3"),

            br(),

            plotOutput("mPlot3", height = "550px"),

            br(),
            br(),

            strong('R session info'),
            verbatimTextOutput("info3.out")
            ),









        tabPanel("Test of Independence (Tabulated data)",

            h2("Test of Independence (Tabulated data)"),

            h4("Two or more than two nominal variables"),

            p('Note: Input values must be separated by tabs. Copy and paste from Excel/Numbers.'),

            p(HTML("<b><div style='background-color:#FADDF2;border:1px solid black;'>Your data needs to have the header (variable names) in the first row. Missing values should be indicated by a period (.) or NA.</div></b>")),

            aceEditor("text4", value="\tNo\tYes\nM\t20\t18\nW\t8\t24", mode="r", theme="cobalt"),

            br(),

            h3("Contingency table"),
            verbatimTextOutput("data4.out"),

            br(),

            h3("Test result"),
            verbatimTextOutput("test4.out"),

            br(),

            h3("Plot"),

            plotOutput("pPlot4"),

            br(),

            plotOutput("mPlot4", height = "550px"),

            br(),
            br(),

            strong('R session info'),
            verbatimTextOutput("info4.out")
            ),










        tabPanel("About",

            strong('Note'),
            p('This web application is developed with',
            a("Shiny.", href="http://www.rstudio.com/shiny/", target="_blank"),
            ''),

            br(),

            strong('List of Packages Used'), br(),
            code('library(shiny)'),br(),
            code('library(shinyAce)'),br(),
            code('library(pwr)'),br(),
            code('library(vcd)'),br(),

            br(),

            strong('Code'),
            p('Source code for this application is based on',
            a('"The handbook of Research in Foreign Language Learning and Teaching" (Takeuchi & Mizumoto, 2012).', href='http://mizumot.com/handbook/', target="_blank")),

            p('The code for this web application is available at',
            a('GitHub.', href='https://github.com/mizumot/chi', target="_blank")),

            p('If you want to run this code on your computer (in a local R session), run the code below:',
            br(),
            code('library(shiny)'),br(),
            code('runGitHub("chi","mizumot")')
            ),

            p('I referred to',
            a("js-STAR", href="http://www.kisnet.or.jp/nappa/software/star/", target="_blank"),
            'for some parts of the codes. I would like to thank the authors of js-STAR, the very fast and excellent online software.'),

            br(),

            strong('Citation in Publications'),
            p('Mizumoto, A. (2015). Langtest (Version 1.0) [Web application]. Retrieved from http://langtest.jp'),

            br(),

            strong('Article'),
            p('Mizumoto, A., & Plonsky, L. (2015).', a("R as a lingua franca: Advantages of using R for quantitative research in applied linguistics.", href='http://applij.oxfordjournals.org/content/early/2015/06/24/applin.amv025.abstract', target="_blank"), em('Applied Linguistics,'), 'Advance online publication. doi:10.1093/applin/amv025'),

            br(),

            strong('Recommended'),
            p('To learn more about R, I suggest this excellent and free e-book (pdf),',
            a("A Guide to Doing Statistics in Second Language Research Using R,", href="http://cw.routledge.com/textbooks/9780805861853/guide-to-R.asp", target="_blank"),
            'written by Dr. Jenifer Larson-Hall.'),

            p('Also, if you are a cool Mac user and want to use R with GUI,',
            a("MacR", href="https://sites.google.com/site/casualmacr/", target="_blank"),
            'is defenitely the way to go!'),

            br(),

            strong('Author'),
            p(a("Atsushi MIZUMOTO,", href="http://mizumot.com", target="_blank"),' Ph.D.',br(),
            'Associate Professor of Applied Linguistics',br(),
            'Faculty of Foreign Language Studies /',br(),
            'Graduate School of Foreign Language Education and Research,',br(),
            'Kansai University, Osaka, Japan'),

            br(),

            a(img(src="http://i.creativecommons.org/p/mark/1.0/80x15.png"), target="_blank", href="http://creativecommons.org/publicdomain/mark/1.0/"),

            p(br())

)
)
)
))
Code license: GPL-3