# R analysis for the manuscript
# Schmid, D., Hesse, C., & Schenk, T. (in review). Blindsight and residual vision: How reliable is the redundant target effect as a
#  tool? Cortex


# using R version 4.2.0 (2022-04-22 ucrt) -- "Vigorous Calisthenics"
# calculate 95% confidence intervals using the Wilson score method

# use binomial package
# https://www.rdocumentation.org/packages/binom/versions/1.1-1.1/topics/binom.confint
# install.packages("binom")   # run once
require(binom)
require(tidyverse)
require(readxl)

# set current directory
# current directory has to contain the following files
# RTP Review.R
# RTP Review Data.xlsx
setwd("XXX") # insert your directory here!

# load data from excel file
Data <- read_xlsx("RTP Review Data.xlsx")

# prevalence of RVCs in patients combining single-case and group results: 21/74; see table 1
X = sum(Data$`RVC-Prevalence N Clean`, na.rm = TRUE)
N = sum(Data$`HVFD Patients N Clean`, na.rm = TRUE)
result <- binom.confint(X, N, methods = "wilson")
print(result)
paste("value = ",round(result$mean*100, digits = 1), "% (95% CI [", round(result$lower*100),"%, ",round(result$upper*100),"%])", sep = "")

# prevalence of RVCs in patients based on single-case results: 15/68
X = sum(Data$`RVC-Prevalence SC N Clean`, na.rm = TRUE)
N = sum(Data$`RVC SC-Test N Clean`, na.rm = TRUE)
result <- binom.confint(X, N, methods = "wilson")
print(result)
paste("value = ",round(result$mean*100, digits = 1), "% (95% CI [", round(result$lower*100),"%, ",round(result$upper*100),"%])", sep = "")

# sensitivity is assessed by the RTE measured for a unilateral-redundant condition where targets are presented within the sighted visual field
# sensitivity in patients based on single-case results: 5/12
X = sum(Data$`RTE Unilateral SC N`, na.rm = TRUE)
N = sum(Data$`Unilateral SC N`, na.rm = TRUE)
result <- binom.confint(X, N, methods = "wilson")
print(result)
paste("value = ",round(result$mean*100, digits = 1), "% (95% CI [", round(result$lower*100),"%, ",round(result$upper*100),"%])", sep = "")

# sensitivity in patients combining single-case and group results: 34/53
X = sum(Data$`RTE Unilateral N`, na.rm = TRUE)
N = sum(Data$`Unilateral Test N`, na.rm = TRUE)
result <- binom.confint(X, N, methods = "wilson")
print(result)
paste("value = ",round(result$mean*100, digits = 1), "% (95% CI [", round(result$lower*100),"%, ",round(result$upper*100),"%])", sep = "")

# sensitivity in healthy participants: 17/22
# derived from Schärli et al., 1999
X = sum(Data[Data$Study == "Schärli et al., 1999", "Healthy Participants RTE-Prevalence"])
N = sum(Data[Data$Study == "Schärli et al., 1999", "Healthy Participants SC N"])
result <- binom.confint(X, N, methods = "wilson")
print(result)
paste("value = ",round(result$mean*100, digits = 1), "% (95% CI [", round(result$lower*100),"%, ",round(result$upper*100),"%])", sep = "")

# specificity in healthy participants: 17/19 do NOT show the effect
# derived from Schmid et al., 2022
X = sum(Data[Data$Study == "Schmid et al., 2022", "Healthy Participants RTE-Absence"])
N = sum(Data[Data$Study == "Schmid et al., 2022", "Healthy Participants SC N"])
result <- binom.confint(X, N, methods = "wilson")
print(result)
paste("value = ",round(result$mean*100, digits = 1), "% (95% CI [", round(result$lower*100),"%, ",round(result$upper*100),"%])", sep = "")

# reliability in patients: 4/11; see Table 3
X = sum(Data$`RTE Reliable N Clean`, na.rm = TRUE)
N = sum(Data$`Repeated Testing N Clean`, na.rm = TRUE)
result <- binom.confint(x = X, n = N, methods = "wilson")
print(result)
paste("value = ",round(result$mean*100, digits = 1), "% (95% CI [", round(result$lower*100),"%, ",round(result$upper*100),"%])", sep = "")

# unreliable results in patients: 7/11; see Table 3
X = sum(Data$`Repeated Testing N Clean`, na.rm = TRUE) - sum(Data$`RTE Reliable N Clean`, na.rm = TRUE)
N = sum(Data$`Repeated Testing N Clean`, na.rm = TRUE)
result <- binom.confint(x = X, n = N, methods = "wilson")
print(result)
paste("value = ",round(result$mean*100, digits = 1), "% (95% CI [", round(result$lower*100),"%, ",round(result$upper*100),"%])", sep = "")


# -------------------------------------------------------------
# create figure for positive and negative predictive values

# load packages
library(tidyverse)
library(ggpattern)

# define folder to save figures
# in this folder, there should also be the PlotSytleApa.R file
setwd("E:/Doris/Study RTP Review/8-Upload/")

# load default plot style
source("PlotStyleApa.R")

# prevalence estimate based on Marzi et al., 1986: 4/20 = 20%
Prevalence = 20 # 20 of 100 patients possess RVCs
Sensitivity = 0.42
Specificity = 0.89
TruePos20 = round(Prevalence*Sensitivity) # RVC cases correctly detected
FalseNeg20 = Prevalence-TruePos20           # RVC cases missed
TrueNeg20 = round((100-Prevalence)*Specificity) # No-RVC cases correctly detected
FalsePos20 = (100-Prevalence) - TrueNeg20

AllPos20 = TruePos20 + FalsePos20
AllNeg20 = TrueNeg20 + FalseNeg20

# prevalence estimate: 60%
Prevalence = 60 # 60 of 100 patients possess RVCs
Sensitivity = 0.42
Specificity = 0.89
TruePos60 = round(Prevalence*Sensitivity) # RVC cases correctly detected
FalseNeg60 = Prevalence-TruePos60           # RVC cases missed
TrueNeg60 = round((100-Prevalence)*Specificity) # No-RVC cases correctly detected
FalsePos60 = (100-Prevalence) - TrueNeg60

AllPos60 = TruePos60 + FalsePos60
AllNeg60 = TrueNeg60 + FalseNeg60

# define data
Values <- data.frame(Prevalence = c("Prevalence 20%","Prevalence 20%", "Prevalence 20%","Prevalence 20%", "Prevalence 60%", "Prevalence 60%", "Prevalence 60%", "Prevalence 60%"),
                     Outcome = c('Positive','Negative','Positive','Negative','Positive','Negative','Positive','Negative'),
                     ConditionNumber = c(AllPos20, AllNeg20, AllPos20, AllNeg20, AllPos60, AllNeg60, AllPos60, AllNeg60),
                     Correct = c('Correct','Correct','False','False','Correct','Correct','False','False'),
                     Count = c(TruePos20, TrueNeg20, FalsePos20, FalseNeg20, TruePos60, TrueNeg60, FalsePos60,FalseNeg60))

Values$Percent <- round(Values$Count/Values$ConditionNumber*100)
Values$PercentStr <- paste(Values$Percent,"%", sep= "")
Values$PercentY <- (Values$ConditionNumber-(Values$ConditionNumber-Values$Count))/2 +(Values$ConditionNumber-Values$Count)
Values$PercentY[Values$Correct == "False"] = Values$Count[Values$Correct == "False"]/2

Values$Outcome <- factor(Values$Outcome, levels = c("Positive","Negative"))

# define figure name for png
figurename_png = paste("RTP_Figure_1.png", sep = "")

# define size of png-figure
# 1 column: 9.14 cm; 1.5 column: 14.09 cm; 2 columns: 19.05 cm
width_cm = 14.09
height_cm = 9.14

print(ggplot(data = Values, aes(x=Outcome, y = Count, fill = Correct))+              
        geom_col(color = "black")+
        scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 25))+
        scale_fill_manual(values = c("white","grey"))+
        geom_text(data = Values, aes(x = Outcome, y = PercentY, label = Count), color = "black", size = 2.6)+
        ylab("Number of Cases")+
        xlab("Test Outcome")+
        facet_grid(cols = vars(Prevalence))+
        PlotStyle+
        labs(fill = NULL))         

# save figure as png
ggsave(filename = figurename_png, width = width_cm, height = height_cm, units = "cm", dpi = 300)