Comparing Paired Samples

library(shiny)
library(psych)
library(car)
library(compute.es)
library(pwr)
library(lattice)
library(latticeExtra)
library(beeswarm)

shinyServer(function(input, output) {
    options(warn = -1)

    # データの取得と前処理
    getData <- reactive({
        if (input$colname) {
            # ヘッダーありの場合
            df <- read.csv(text = input$text, sep = "\t", na.strings = c("", "NA", ".", "#N/A", "N/A"))
            # 完全なケースのみを抽出（リストワイズ削除）
            df_complete <- df[complete.cases(df[, 1:2]), ]
            x <- df_complete[, 1] # 1列目のデータ
            y <- df_complete[, 2] # 2列目のデータ
            var_names <- colnames(df)[1:2] # 変数名を取得
        } else {
            # ヘッダーなしの場合
            df <- read.csv(text = input$text, header = FALSE, sep = "\t", na.strings = c("", "NA", ".", "#N/A", "N/A"))
            # 完全なケースのみを抽出（リストワイズ削除）
            df_complete <- df[complete.cases(df[, 1:2]), ]
            x <- df_complete[, 1] # 1列目のデータ
            y <- df_complete[, 2] # 2列目のデータ
            var_names <- c("Data 1", "Data 2") # デフォルト変数名
        }

        # 欠損値の情報を追加
        total_rows <- nrow(df)
        complete_rows <- nrow(df_complete)
        missing_rows <- total_rows - complete_rows

        return(list(
            x = x,
            y = y,
            df = df_complete,
            var_names = var_names,
            total_rows = total_rows,
            complete_rows = complete_rows,
            missing_rows = missing_rows
        ))
    })

    bs <- reactive({
        data <- getData()
        x <- data$x
        y <- data$y
        var_names <- data$var_names

        # 欠損値の情報
        total_rows <- data$total_rows
        complete_rows <- data$complete_rows
        missing_rows <- data$missing_rows

        result1 <- describe(x)[2:13]
        result2 <- describe(y)[2:13]
        row.names(result1) <- ""
        row.names(result2) <- ""

        correl <- cor(x, y)

        result <- list()
        result[[var_names[1]]] <- result1
        result[[var_names[2]]] <- result2
        result[["Correlation"]] <- correl

        # 欠損値の情報を追加
        if (missing_rows > 0) {
            result[["Missing Values"]] <- c(
                paste("Total rows:", total_rows),
                paste("Complete cases:", complete_rows),
                paste("Missing cases (list-wise deleted):", missing_rows)
            )
        }

        return(result)
    })


    makedistPlot <- function() {
        data <- getData()
        x <- data$x
        y <- data$y
        var_names <- data$var_names

        simple.bincount <- function(x, breaks) {
            nx <- length(x)
            nbreaks <- length(breaks)
            counts <- integer(nbreaks - 1)
            for (i in 1:nx) {
                lo <- 1
                hi <- nbreaks
                if (breaks[lo] <= x[i] && x[i] <= breaks[hi]) {
                    while (hi - lo >= 2) {
                        new <- (hi + lo) %/% 2
                        if (x[i] > breaks[new]) {
                            lo <- new
                        } else {
                            hi <- new
                        }
                    }
                    counts[lo] <- counts[lo] + 1
                }
            }
            return(counts)
        }

        nclass.x <- nclass.FD(x)
        breaks.x <- pretty(x, nclass.x)
        counts.x <- simple.bincount(x, breaks.x)
        counts.max.x <- max(counts.x)

        nclass.y <- nclass.FD(y)
        breaks.y <- pretty(y, nclass.y)
        counts.y <- simple.bincount(y, breaks.y)
        counts.max.y <- max(counts.y)

        counts.max <- max(c(counts.max.x, counts.max.y))


        xy.min <- min(c(x, y))
        xy.min <- xy.min - xy.min * 0.1
        xy.max <- max(c(x, y))
        xy.max <- xy.max + xy.max * 0.1

        p1 <- hist(x, xlim = c(xy.min, xy.max), ylim = c(0, counts.max * 1.3))
        p2 <- hist(y, xlim = c(xy.min, xy.max), ylim = c(0, counts.max * 1.3))

        plot(p1,
            las = 1, xlab = paste(var_names[1], "is expressed in blue;", var_names[2], "in red. Vertical lines show the mean."),
            main = "", col = rgb(0, 0, 1, 1 / 4), xlim = c(xy.min, xy.max), ylim = c(0, counts.max * 1.3)
        )
        plot(p2, las = 1, xlab = "", main = "", col = rgb(1, 0, 0, 1 / 4), xlim = c(xy.min, xy.max), ylim = c(0, counts.max * 1.3), add = TRUE)

        abline(v = mean(x), col = "blue", lwd = 2)
        abline(v = mean(y), col = "red", lwd = 2)
    }

    output$distPlot <- renderPlot({
        print(makedistPlot())
    })



    makeboxPlot <- function() {
        data <- getData()
        x <- data$x
        y <- data$y
        var_names <- data$var_names

        score <- c(x, y)
        group <- factor(c(rep(var_names[1], length(x)), rep(var_names[2], length(y))))

        boxplot(score ~ group, las = 1, xlab = "Means and +/-1 SDs are displayed in red.")

        beeswarm(score ~ group, col = 4, pch = 16, add = TRUE)

        points(1.2, mean(x), pch = 18, col = "red", cex = 2)
        arrows(1.2, mean(x), 1.2, mean(x) + sd(x), length = 0.1, angle = 45, col = "red")
        arrows(1.2, mean(x), 1.2, mean(x) - sd(x), length = 0.1, angle = 45, col = "red")

        points(2.2, mean(y), pch = 18, col = "red", cex = 2)
        arrows(2.2, mean(y), 2.2, mean(y) + sd(y), length = 0.1, angle = 45, col = "red")
        arrows(2.2, mean(y), 2.2, mean(y) - sd(y), length = 0.1, angle = 45, col = "red")
    }

    output$boxPlot <- renderPlot({
        print(makeboxPlot())
    })



    makeindvPlot <- function() {
        data <- getData()
        x <- data$x
        y <- data$y
        var_names <- data$var_names

        xy <- c(x, y)
        data.point <- factor(c(rep(var_names[1], length(x)), rep(var_names[2], length(y))))
        dat <- data.frame(xy, data.point)
        dat$indiv <- factor(c(rep(1:length(x)), rep(1:length(y))))
        each <- xyplot(xy ~ data.point,
            group = indiv,
            type = c("l"), data = dat, xlab = "", ylab = ""
        )
        a <- mean(x)
        b <- mean(y)
        value <- c(a, b)
        data.point2 <- factor(c(rep(var_names[1], 1), rep(var_names[2], 1)))
        dat2 <- data.frame(value, data.point2)
        all <- xyplot(value ~ data.point2,
            col = "black",
            lwd = 5, type = c("l"), data = dat, xlab = "", ylab = ""
        )
        indv <- each + as.layer(all, axes = NULL)
        print(indv)
    }

    output$indvPlot <- renderPlot({
        print(makeindvPlot())
    })



    makecorrelPlot <- function() {
        data <- getData()
        x <- data$x
        y <- data$y
        var_names <- data$var_names

        xy.min <- min(c(x, y))
        xy.min <- xy.min - xy.min * 0.1
        xy.max <- max(c(x, y))
        xy.max <- xy.max + xy.max * 0.1

        plot(x, y,
            las = 1, pch = 16, xlab = var_names[1], ylab = var_names[2], main = "",
            xlim = c(xy.min, xy.max), ylim = c(xy.min, xy.max)
        )
        lines(par()$usr[1:2], par()$usr[3:4], lty = 3)
    }

    output$correlPlot <- renderPlot({
        print(makecorrelPlot())
    })



    testnorm <- reactive({
        data <- getData()
        x <- data$x
        y <- data$y
        var_names <- data$var_names

        data.1ks <- ks.test(scale(x), "pnorm")
        data.1sh <- shapiro.test(x)

        data.2ks <- ks.test(scale(y), "pnorm")
        data.2sh <- shapiro.test(y)

        result <- list()
        result[[var_names[1]]] <- data.1ks
        result[[paste0(var_names[1], " (n < 50 w/ Outliers)")]] <- data.1sh
        result[[var_names[2]]] <- data.2ks
        result[[paste0(var_names[2], " (n < 50 w/ Outliers)")]] <- data.2sh

        return(result)
    })



    makedistdiffPlot <- function() {
        data <- getData()
        x <- data$x
        y <- data$y

        diff <- x - y

        simple.bincount <- function(x, breaks) {
            nx <- length(diff)
            nbreaks <- length(breaks)
            counts <- integer(nbreaks - 1)
            for (i in 1:nx) {
                lo <- 1
                hi <- nbreaks
                if (breaks[lo] <= x[i] && x[i] <= breaks[hi]) {
                    while (hi - lo >= 2) {
                        new <- (hi + lo) %/% 2
                        if (x[i] > breaks[new]) {
                            lo <- new
                        } else {
                            hi <- new
                        }
                    }
                    counts[lo] <- counts[lo] + 1
                }
            }
            return(counts)
        }

        nclass <- nclass.FD(diff)
        breaks <- pretty(diff, nclass)
        counts <- simple.bincount(diff, breaks)
        counts.max <- max(counts)

        h <- hist(diff,
            las = 1, breaks = "FD", xlab = "Red vertical line shows the mean of the differences.",
            ylim = c(0, counts.max * 1.2), main = "", col = "cyan"
        )
        rug(diff)
        abline(v = mean(diff, na.rm = T), col = "red", lwd = 2)
        xfit <- seq(min(diff, na.rm = T), max(diff, na.rm = T))
        yfit <- dnorm(xfit, mean = mean(diff, na.rm = T), sd = sd(diff, na.rm = T))
        yfit <- yfit * diff(h$mids[1:2]) * length(diff)
        lines(xfit, yfit, col = "blue", lwd = 2)
    }

    output$distdiffPlot <- renderPlot({
        print(makedistdiffPlot())
    })



    diffnorm <- reactive({
        data <- getData()
        x <- data$x
        y <- data$y

        diff <- x - y

        diff.1ks <- ks.test(scale(diff), "pnorm")
        diff.1sh <- shapiro.test(diff)

        return(list(Diff.Score = diff.1ks, "Diff.Score (n < 50 w/ Outliers)" = diff.1sh))
    })



    t <- reactive({
        data <- getData()
        x <- data$x
        y <- data$y

        t.test(x, y, paired = TRUE)
    })



    es <- reactive({
        data <- getData()
        x <- data$x
        y <- data$y

        m1 <- mean(x)
        sd1 <- sd(x)
        n1 <- length(x)

        m2 <- mean(y)
        sd2 <- sd(y)
        n2 <- length(y)

        diff <- x - y

        pooledSD <- sd(diff)
        sediff <- pooledSD / sqrt(length(x))
        paired.d <- abs((m1 - m2) / pooledSD)
        dCI <- d.ci(paired.d, n1 = n1, alpha = .05)

        d <- abs(mean(diff) / (sd(diff) / sqrt(2 * (1 - cor(x, y)))))
        var.d <- ((1 / n1) + ((d^2) / (2 * n1))) * (2 * (1 - cor(x, y)))

        df <- n1 - 1
        j <- 1 - (3 / (4 * df - 1))
        g <- j * d
        var.g <- j^2 * var.d

        # a <- ((n1 + n2)^2)/(n1 * n2)
        # r <- d/sqrt((d^2) + a)
        # var.r <- (a^2 * var.d)/(d^2 + a)^3

        alpha <- 0.05
        crit <- qt(alpha / 2, df, lower.tail = FALSE)

        lower.d <- d - crit * sqrt(var.d)
        upper.d <- d + crit * sqrt(var.d)

        # lower.d2 <- d - 1.96 * sqrt(var.d)
        # upper.d2 <- d + 1.96 * sqrt(var.d)

        lower.g <- g - crit * sqrt(var.g)
        upper.g <- g + crit * sqrt(var.g)

        # lower.g2 <- d - 1.96 * sqrt(var.d)
        # upper.g2 <- d + 1.96 * sqrt(var.d)

        cat("=======================================================", "\n")
        cat(
            " paired.t / sqrt(n):", "\n",
            "(Based on Larkens, 2013, p. 4)", "\n",
            "\n",
            "  d [95% CI] =", dCI[2], "[", dCI[1], ",", dCI[3], "]", "\n",
            "\n"
        )



        cat("=======================================================", "\n")
        cat(
            " Mean difference / within-groups SD:", "\n",
            "(Based on Borenstein et al., 2021, p. 29) ", "\n",
            "\n",
            "  d [95% CI] =", d, "[", lower.d, ",", upper.d, "]", "\n",
            "   var(d) =", var.d, "\n",
            "\n",
            "  g [95% CI] =", g, "[", lower.g, ",", upper.g, "]", "\n",
            "   var(g) =", var.g, "\n",
            "\n"
        )



        delta <- abs(mean(diff) / sd(x))

        cat("=======================================================", "\n")
        cat(
            " Mean difference / SD of Data 1 (e.g., pretest):", "\n",
            "\n",
            "  Delta =", delta, "\n",
            "\n"
        )



        result.dependent <- t.test(x, y, paired = TRUE)
        paired.t <- result.dependent$statistic
        r2 <- sqrt(paired.t^2 / (paired.t^2 + df))
        r2 <- r2[[1]]

        cat("=======================================================", "\n")
        cat(
            " r --- sqrt(paired.t^2/(paired.t^2+df)):", "\n",
            "\n",
            "  r =", r2, "\n",
            "\n"
        )
    })



    wsr <- reactive({
        data <- getData()
        x <- data$x
        y <- data$y

        result <- wilcox.test(x, y, paired = TRUE, correct = FALSE)

        pval <- result$p.value
        z <- qnorm(1 - (pval / 2))

        r1 <- z / sqrt(length(x * 2))

        r2 <- z / sqrt(length(x) - sum((y - x == 0)))

        print(result)

        cat(
            " z-value =", z, "\n",
            "effect size r =", r1, "\n",
            "effect size r (without considering ties) =", r2, "\n"
        )
    })



    power <- reactive({
        data <- getData()
        x <- data$x
        y <- data$y

        diff <- x - y
        d <- mean(diff) / sd(diff)
        n <- length(x)

        posthoc <- power.t.test(n = n, delta = d, sig.level = .05, type = "paired", strict = TRUE)$power

        future <- ceiling(power.t.test(
            power = 0.8, delta = d, sig.level = 0.05,
            type = "paired", strict = TRUE, alternative = "two.sided"
        )$n)

        cat(
            " Post hoc (observed) power =", round(posthoc, 3), "\n",
            "\n",
            "  Note: According to Cumming (2012), post hoc power is 'illegitimate'", "\n",
            "        and we should NEVER calculate or report it.", "\n",
            "\n",
            "\n",
            "Sample size needed for future experiment:", " n =", future, "\n",
            "  Power = 0.8, sig.level = 0.05, alternative = two.sided, d =", round(d, 2), "\n",
            "\n",
            "  Note: This is true only PROVIDED that the population effect size is", "\n",
            "        equal to the observed sample effect size (i.e., it is unrealistic).", "\n",
            "\n",
            "\n",
            "POWER ANALYSIS SHOULD BE CONDUCTED PRIOR TO THE EXPERIMENT!", "\n"
        )
    })



    info <- 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$info.out <- renderPrint({
        info()
    })



    output$textarea.out <- renderPrint({
        bs()
    })

    output$testnorm.out <- renderPrint({
        testnorm()
    })

    output$diffnorm.out <- renderPrint({
        diffnorm()
    })

    output$t.out <- renderPrint({
        t()
    })

    output$es.out <- renderPrint({
        es()
    })

    output$wsr.out <- renderPrint({
        wsr()
    })

    output$power.out <- renderPrint({
        power()
    })
})