## ------------------------------------------------------------------------
library(sensR)
discrim(15, 20, method = "tetrad", conf.level = 0.90)

## ------------------------------------------------------------------------
?discrim

## ------------------------------------------------------------------------
library(binom)
binom.confint(15, n = 20, conf.level = 0.90, methods = "exact")

## ------------------------------------------------------------------------
pc2pd(c(0.75, 0.5444176, 0.8959192), Pguess = 1/3)

## ------------------------------------------------------------------------
psyinv(c(0.75, 0.5444176, 0.8959192), method = "tetrad")

## ------------------------------------------------------------------------
rescale(pc = c(0.75, 0.5444176, 0.8959192), method = "tetrad")
rescale(d.prime =  c(1.889770, 1.176030, 2.604497), method = "tetrad")
rescale(pd = c(0.6250000, 0.3166264, 0.8438788), method = "tetrad")

## ------------------------------------------------------------------------
findcr(20, alpha = 0.05, p0 = 1/3)

## ------------------------------------------------------------------------
findcr(20, alpha = 0.05, p0 = 1/3, pd0 = 0.50, test = "similarity")

## ------------------------------------------------------------------------
findcr(20, alpha = 0.05, p0 = 1/3, pd0 = pc2pd(psyfun(1.5, "tetrad"), Pguess=1/3), 
       test = "similarity")

## ------------------------------------------------------------------------
pc2pd(psyfun(1.5, "tetrad"), Pguess=1/3)

## ------------------------------------------------------------------------
discrim(28, 60, d.prime0 = 1.3, method = "tetrad", conf.level = 0.90, 
        test = "similarity")

## ------------------------------------------------------------------------
plot(discrim(15, 20, method = "tetrad"), cex.main=0.8)

## ------------------------------------------------------------------------
CIres <- binom.confint(c(15, 7, 12, 6), 20, conf.level = 0.90, methods = "exact")

## ------------------------------------------------------------------------
pc2pd(as.matrix(CIres[,4:6]), Pguess = 1/10)

## ------------------------------------------------------------------------
binom.test(6, 20, p=1/10, alternative = "greater")

## ------------------------------------------------------------------------
discrimPwr(pdA = 0.5, sample.size = 20, alpha = 0.05, pGuess = 1/3)

## ------------------------------------------------------------------------
d.primePwr(d.primeA = 1.5, sample.size = 20, alpha = 0.05, method = "tetrad")

## ------------------------------------------------------------------------
d.primePwr(d.primeA = 1.5, sample.size = 20, alpha = 0.05, method = "duotrio")
d.primePwr(d.primeA = 1.5, sample.size = 20, alpha = 0.05, method = "triangle")
d.primePwr(d.primeA = 1.5, sample.size = 20, alpha = 0.05, method = "tetrad")
d.primePwr(d.primeA = 1.5, sample.size = 20, alpha = 0.05, method = "twoAFC")
d.primePwr(d.primeA = 1.5, sample.size = 20, alpha = 0.05, method = "threeAFC")

## ------------------------------------------------------------------------
# Instead of these simple five calls the same is done slightly different to 
# make a nicer one line and rounded printout
protocols <- c("duotrio", "triangle", "tetrad",  "twoAFC", "threeAFC")
powers <- rep(0, 5)
for (i in 1:5) powers[i] <- d.primePwr(d.primeA = 1.5, sample.size = 20, 
                                       alpha = 0.05, method = protocols[i])
names(powers) <- protocols
round(powers, 4)

## ------------------------------------------------------------------------
discrimPwr(pdA = pc2pd(psyfun(1.5, "duotrio", double = TRUE), Pguess = 1/4),
            sample.size = 20, alpha = 0.05, pGuess = 1/4)
discrimPwr(pdA = pc2pd(psyfun(1.5, "triangle", double = TRUE), Pguess = 1/9),
            sample.size = 20, alpha = 0.05, pGuess = 1/9)
discrimPwr(pdA = pc2pd(psyfun(1.5, "twoAFC", double = TRUE), Pguess = 1/4),
            sample.size = 20, alpha = 0.05, pGuess = 1/4)
discrimPwr(pdA = pc2pd(psyfun(1.5, "threeAFC", double = TRUE), Pguess = 1/9),
            sample.size = 20, alpha = 0.05, pGuess = 1/9)

## ------------------------------------------------------------------------
# Instead of these simple four calls the same is done slightly different to 
# make a nicer one line and rounded printout
protocols <- c("duotrio", "triangle", "twoAFC", "threeAFC")
pGuess <- c(1/4, 1/9, 1/4, 1/9)
powers <- rep(0, 4)
for (i in 1:4) powers[i] <- discrimPwr(pdA = pc2pd(psyfun(1.5, protocols[i], double = TRUE), 
        Pguess=pGuess[i]), sample.size = 20, alpha = 0.05, pGuess = pGuess[i])
names(powers) <- paste0("double_",protocols)
round(powers, 4)

## ------------------------------------------------------------------------
## R-code to produce Figure 15.3 (not shown in the chapter)
n <- 20
dAs <- seq(0, 3, 0.05)
powers <- dAs
doub_tri_powers <- powers

tetrad_powers <- powers
duotrio_powers  <- powers
doub_duotrio_powers  <- powers

twoAFC_powers  <- powers
doub_2AFC_powers  <- powers
threeAFC_powers  <- powers
doub_3AFC_powers  <- powers


for (i in 1:length(dAs)) powers[i] <-
  d.primePwr(dAs[i], method="triangle", sample.size = n)
for (i in 1:length(dAs)) tetrad_powers[i] <-
  d.primePwr(dAs[i], method="tetrad", sample.size = n)
for (i in 1:length(dAs)) duotrio_powers[i] <-
  d.primePwr(dAs[i], method="duotrio", sample.size = n)
for (i in 1:length(dAs)) twoAFC_powers[i] <-
  d.primePwr(dAs[i], method="twoAFC", sample.size = n)
for (i in 1:length(dAs)) threeAFC_powers[i] <-
  d.primePwr(dAs[i], method="threeAFC", sample.size = n)

for (i in 1:length(dAs)) doub_3AFC_powers[i] <-
  discrimPwr(pc2pd(psyfun(dAs[i], method = "threeAFC", double = TRUE), 1/9),
             pGuess = 1/9, sample.size = n)
for (i in 1:length(dAs)) doub_tri_powers[i] <-
  discrimPwr(pc2pd(psyfun(dAs[i], method = "triangle", double = TRUE), 1/9),
             pGuess = 1/9, sample.size = n)

for (i in 1:length(dAs)) doub_duotrio_powers[i] <-
  discrimPwr(pc2pd(psyfun(dAs[i], method = "duotrio", double = TRUE), 1/4),
             pGuess = 1/4, sample.size = n)
for (i in 1:length(dAs)) doub_2AFC_powers[i] <-
  discrimPwr(pc2pd(psyfun(dAs[i], method = "twoAFC", double = TRUE), 1/4),
             pGuess = 1/4, sample.size = n)


plot(dAs, powers, type="n", xlab = "d-prime", ylab = "Power", main = "n=20")

lines(dAs, powers, col="black")
lines(dAs, doub_tri_powers, col="black", lty = 3)

lines(dAs, duotrio_powers, col="green")
lines(dAs, doub_duotrio_powers, col="green", lty = 3)

lines(dAs, twoAFC_powers, col="blue")
lines(dAs, doub_2AFC_powers, col="blue", lty = 3)

lines(dAs, threeAFC_powers, col="red")
lines(dAs, doub_3AFC_powers, col="red", lty = 3)

lines(dAs, tetrad_powers, col="grey")

i1 <- 5:8
i2 <- c(1:4,9)
legends <- c("Triangle", "Doub_triangle", "Duotrio", "Doub_duotrio",
                               "twoAFC", "Doub_2AFC","threeAFC", "Doub_3AFC","Tetrad")
cols <- c("black", "black", "green", "green", "blue", "blue", "red", "red", "grey")
ltys <- c(1, 3, 1, 3, 1, 3, 1, 3, 1)
legend("topleft", legend=legends[i1], text.col = cols[i1], lty = ltys[i1], col = cols[i1])
legend("bottomright", legend=legends[i2], text.col = cols[i2], lty = ltys[i2], col = cols[i2])

## ------------------------------------------------------------------------
d.primePwr(d.primeA = 0, d.prime0 = 1, sample.size = 20, alpha = 0.05, 
           method = "tetrad", test = "similarity")
d.primePwr(d.primeA = 0, d.prime0 = 1, sample.size = 100, alpha = 0.05, 
           method = "tetrad", test = "similarity")

## ------------------------------------------------------------------------
d.primePwr(d.primeA = 0.5, d.prime0 = 1, sample.size = 100, alpha = 0.05, 
           method = "tetrad", test = "similarity")

## ------------------------------------------------------------------------
d.primeSS(d.primeA = 0, d.prime0 = 1, target.power = 0.9, alpha = 0.05, 
          method = "tetrad", test = "similarity")
d.primeSS(d.primeA = 0.5, d.prime0 = 1, target.power = 0.9, alpha = 0.05, 
          method = "tetrad", test = "similarity")

## ------------------------------------------------------------------------
d.primeSS(d.primeA = 1, target.power = 0.9, method = "tetrad")
discrimSS(pdA = 0.2407, target.power = 0.9, pGuess = 1/3)
rescale(d.prime=1, method="tetrad")

## ------------------------------------------------------------------------
d.primeSS(d.primeA = 1, target.power = 0.9, alpha = 0.20, method = "tetrad")
d.primeSS(d.primeA = 1, target.power = 0.9, alpha = 0.20, method = "tetrad",
          statistic = "stable.exact")

## ------------------------------------------------------------------------
1-d.primePwr(1, sample.size = 49, method = "tetrad", alpha=0.2)

## ------------------------------------------------------------------------
findcr(49, alpha = 0.2, p0 = 1/3 )

## ------------------------------------------------------------------------
pc2pd(psyfun(1, method = "tetrad"), 1/3)
findcr(49, alpha = 0.1, p0 = 1/3, test = "similarity", pd0 = 0.2407126)

## ------------------------------------------------------------------------
discrim(16, 49, conf.level = 0.80, method = "tetrad")

## ------------------------------------------------------------------------
discrim(16, 49, conf.level = 0.80, method = "tetrad", 
        test="similarity", pd0 = 0.2407126)

## ------------------------------------------------------------------------
findcr(100, alpha = 0.2, p0 = 1/3 )
findcr(100, alpha = 0.1, p0 = 1/3, test = "similarity", pd0 = 0.2407126)

## ------------------------------------------------------------------------
discrim(42, 100, conf.level = 0.80, method = "tetrad", 
        test="similarity", pd0 = 0.2407126)

## ------------------------------------------------------------------------
d.primePwr(d.prime0=1, sample.size = 100, method = "tetrad", alpha=0.1, 
           test="similarity", d.primeA = 0)

## ------------------------------------------------------------------------
discrimSim(15, 10, d.prime = 1, method = "tetrad", sd.indiv = 1)

## ------------------------------------------------------------------------
## R-code to produce Figure 15.4 (not shown in the chapter)
## Set the situation: (make your own choices)
myN=15 #Number of individuals
myK=10 #Number of replications
myd.prime=1 # Level of product difference

par(mfrow = c(2, 2))
par(mar = c(2, 2, 2, 2))

mysd.indiv=0 # individual heterogeneity
# Simulate and plot some data:
rs=discrimSim(myN, replicates = myK, d.prime = myd.prime,method = "tetrad",sd.indiv=mysd.indiv)
barplot(rs[order(rs)],col="blue",ylim=c(0,myK), names.arg=1:myN,cex.names=0.8, main=paste("Heterogeneity:", mysd.indiv))
abline(h=mean(rs))

mysd.indiv=0.5 # individual heterogeneity
# Simulate and plot some data:
rs=discrimSim(myN, replicates = myK, d.prime = myd.prime,method = "tetrad",sd.indiv=mysd.indiv)
barplot(rs[order(rs)],col="blue",ylim=c(0,myK), names.arg=1:myN,cex.names=0.8, main=paste("Heterogeneity:", mysd.indiv))
abline(h=mean(rs))

mysd.indiv=1 # individual heterogeneity
# Simulate and plot some data:
rs=discrimSim(myN, replicates = myK, d.prime = myd.prime,method = "tetrad",sd.indiv=mysd.indiv)
barplot(rs[order(rs)],col="blue",ylim=c(0,myK), names.arg=1:myN,cex.names=0.8, main=paste("Heterogeneity:", mysd.indiv))
abline(h=mean(rs))

mysd.indiv=2 # individual heterogeneity
# Simulate and plot some data:
rs=discrimSim(myN, replicates = myK, d.prime = myd.prime,method = "tetrad",sd.indiv=mysd.indiv)
barplot(rs[order(rs)],col="blue",ylim=c(0,myK), names.arg=1:myN,cex.names=0.8, main=paste("Heterogeneity:", mysd.indiv))
abline(h=mean(rs))
par(mfrow=c(1,1))

## ------------------------------------------------------------------------
#  Replicated data from the Nęs et al book:
x <- c(2, 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 5, 6, 10, 11)
X2 <- cbind(x, 12)

# Analyzed as tetrad data
summary(betabin(X2, method = "tetrad"),  level = 0.9)

## ------------------------------------------------------------------------
1-pchisq(15.492, 2)

## ------------------------------------------------------------------------
(1-pchisq(15.492, 1.5))/2

## ------------------------------------------------------------------------
discrim(sum(x), 180, method = "tetrad", conf.level = 0.9)

## ------------------------------------------------------------------------
qchisq(0.9, 1.5)

## ------------------------------------------------------------------------
## Generate the data:
data <- expand.grid(conc = 1:4, gender = c("Males", "Females"))
data$correct <- c(9, 11, 13, 14, 13, 14, 16, 18)
data$total <- rep(20, 8)
data$concGrp <- factor(data$conc)

## View data:
data

## ------------------------------------------------------------------------
## Fit Initial  model:
glm0 <- glm(cbind(correct, total - correct) ~ gender:concGrp - 1, data,
             family = tetrad)
summary(glm0)

## ------------------------------------------------------------------------
## Fit Inital  model again:
glm0 <- glm(cbind(correct, total - correct) ~ gender*concGrp, data,
             family = tetrad)
anova(glm0, test = "Chisq")

## ------------------------------------------------------------------------
## Fit final model:
glm1 <- glm(cbind(correct, total - correct) ~ gender + conc, data,
             family = tetrad)
## Compare with inital model
anova(glm1, glm0, test = "Chisq")

## ------------------------------------------------------------------------
summary(glm1)

## ------------------------------------------------------------------------
data <- rbind(c(8, 17),
              c(1, 24))
colnames(data) <- c("Yes-responses", "No-responses")
row.names(data) <- c("Yes-samples", "No-samples")
data

## ------------------------------------------------------------------------
H <- 8/(8+17)
FA <- 1/(1+24)
zH <- qnorm(H)
zFA <- qnorm(FA)
## d-prime:
zH - zFA  # d'

## ------------------------------------------------------------------------
data <- rbind(c(8, 17),
              c(1, 24))
SDT(data)

## ------------------------------------------------------------------------
(m1 <- AnotA(8, 25, 1, 25))
confint(m1)[2,]

par(mfrow = c(1, 2))
ROC(m1)
plot(m1, main = "")

## ------------------------------------------------------------------------
# The standard uncorrected pearson test:
chisq.test(matrix(c(21, 21, 9, 33), ncol=2), correct=F)
# The standard corrected pearson test:
chisq.test(matrix(c(21, 21, 9, 33), ncol=2))

## ------------------------------------------------------------------------
qchisq(0.9, 1)

## ------------------------------------------------------------------------
sdres <- samediff(21, 21, 9, 33)
summary(sdres)

# Sensitivity, AUC value with CI:
AUC(sdres$coef[2], sdres$se[2])

## ------------------------------------------------------------------------
pnorm(1.9841/sqrt(2))

## ------------------------------------------------------------------------
samediffPwr(n = 100, tau = 1, delta = 2.5, Ns = 10, Nd = 10)

## ------------------------------------------------------------------------
## Accumulated data:
H <- c(24, 32, 36, 40, 40, 40)/40
FA <-c(1,  2,  6,  8, 16, 36)/36

## The empirical ROC curve:
plot(c(0, FA, 1) , c(0, H, 1))
lines(c(0, FA, 1) , c(0, H, 1))
 
## ------------------------------------------------------------------------
## All the code needed to make Figure 15.6
## Accumulated data:
H <- c(0, 24, 32, 36, 40, 40, 40, 40)/40
FA <-c(0,  1,  2,  6,  8, 16, 36, 36)/36

## The empirical ROC curve:
plot(FA, H)
lines(FA, H)
polygon( c(FA, 1, 0), c(H, 0, 0), col = "grey")
text(0.5, 0.5, "R-index = AUC = 0.953", cex=2)
polygon( c(FA[1:5], 0, 0), c(H[1:5], 1, 0), col = "yellow")
points(FA, H, pch = 16)
lines((0:10)/10, (0:10)/10, lty = 3, lwd=0.5)

## ------------------------------------------------------------------------
## Constructing the two-sample ordinal data from the frequencies
x <- rep(1:6, c(24,8,4,4,0,0))
y <- rep(1:6, c(1,1,4,2,8,20))

## Using inbuilt Wilcoxon function:
(U <- wilcox.test(y, x, correct=F))

## Finding the R-index from that, U/(n1*n2)
(Rindex <- U$statistic/(length(x)*length(y)))

## ------------------------------------------------------------------------
library(ordinal)
# Making the data frame:
mydata <- cbind(c(x, y), c(rep("Same", 40), rep("Diff", 36)))
mydata <- as.data.frame(mydata)
colnames(mydata) <-  c("Score", "Product")                       
mydata$Score <- factor(mydata$Score)

# Analysing by ordinal regression model (cumulative link model, clm):
ordinal_res <- clm(Score ~ Product, data = mydata, link = "probit")
summary(ordinal_res)

# Finding the AUC from the (apparent) d-prime with CI:
AUC(-ordinal_res$coefficients[6], 0.3309)

## ------------------------------------------------------------------------
# Analysing as a Degree-of-difference, DOD data:
dodresults <- dod(mydata$Score[mydata$Product=="Same"], mydata$Score[mydata$Product=="Diff"])
print(dodresults)

## The AUC
AUC(dodresults$d.prime, dodresults$coefficients[2])


## ------------------------------------------------------------------------
dodPwr(d.primeA=1, d.prime0=0, ncat=6, sample.size=100, nsim=1000,
        alpha=.05, method.tau="LR.max", statistic="likelihood")

## ------------------------------------------------------------------------
dodPwr(d.primeA=0, d.prime0=1, ncat=6, sample.size=100, nsim=1000,
        alpha=.05, method.tau="LR.max", statistic="Pearson",
        alternative="similarity")

## ------------------------------------------------------------------------
## Reading the data from a file with same structure as table 1, Chapter 11 
DFCdata1 <- read.table("DFCdata1.txt", header = TRUE, sep =";")
t.test(DFCdata1$Dif)

## ------------------------------------------------------------------------
library(lmerTest)

## Making a "long version"" of the same data 
library(tidyr)
DFCdata1_long <- gather(DFCdata1[,1:3], value = score, 
                        key = Product, Test, Control)
DFCdata1_long$Assessor <- factor(DFCdata1_long$Assessor)
DFCdata1_long$Product <- factor(DFCdata1_long$Product)

## Analysing by mixed model, random assessor effect
lmer1 <- lmer(score ~ Product + (1|Assessor), data = DFCdata1_long)
anova(lmer1)

## ------------------------------------------------------------------------
library(ordinal)
DFCdata1_long$scoreOrd <- ordered(DFCdata1_long$score) 

clmm1 <- clmm(scoreOrd ~ Product + (1|Assessor), data = DFCdata1_long, 
              link = "probit")
summary(clmm1)

## ------------------------------------------------------------------------
clmm2 <- clmm(scoreOrd ~ Product + (1|Assessor), data = DFCdata1_long, 
              link = "probit", threshold = "equidistant")
summary(clmm2)

## ------------------------------------------------------------------------
1-pchisq(2.18, 1)

## ------------------------------------------------------------------------
Controldata <- DFCdata1_long$scoreOrd[DFCdata1_long$Product=="Control"]
Testdata <- DFCdata1_long$scoreOrd[DFCdata1_long$Product=="Test"]
dod(Controldata, Testdata)

## ------------------------------------------------------------------------
## Reading the data from a file with same structure as table 5, Chapter 11
DFCdata2 <- read.table("DFCdata2.txt", header = TRUE, sep =";")
DFCdata2_long <- gather(DFCdata2[,1:4], value = score,
                        key = Product, C, TestA, TestB)
DFCdata2_long$Assessor <- factor(DFCdata2_long$Assessor)
DFCdata2_long$Product <- factor(DFCdata2_long$Product)

## Analysing by mixed model, random assessor effect
lmer2 <- lmer(score ~ Product + (1|Assessor), data = DFCdata2_long)
anova(lmer2)
difflsmeans(lmer2)

## The ordinal analysis
DFCdata2_long$scoreOrd <- factor(DFCdata2_long$score, ordered = TRUE) 

clmm1 <- clmm(scoreOrd ~ Product + (1|Assessor), data = DFCdata2_long, 
              link = "probit")
summary(clmm1)

clmm2 <- clmm(scoreOrd ~ Product + (1|Assessor), data = DFCdata2_long, 
              link = "probit",
              threshold = "equidistant")
summary(clmm2)

## ------------------------------------------------------------------------
1-pchisq(37.68, 1)

## ------------------------------------------------------------------------
summary(lmer2)$sigma

## ------------------------------------------------------------------------
## Reading the data from a file with same structure as table 10, Chapter 12 
Rankdata1 <- read.table("Rankdata1.txt", header = TRUE, sep =";")
Rankdata1_long <- gather(Rankdata1, value = score, 
                         key = Product, Sample1, Sample2, Sample3, Sample4)
Rankdata1_long$Assessor <- factor(Rankdata1_long$Assessor)
Rankdata1_long$Product <- factor(Rankdata1_long$Product)

## Analysing by Friedman's test
(fried1 <- friedman.test(score ~ Product | Assessor, data = Rankdata1_long))

## Doing the post hoc comparison
library(PMCMR)
with(Rankdata1_long, posthoc.friedman.conover.test(score, Product, Assessor, 
                                                   p.adjust="holm"))

## ------------------------------------------------------------------------
findcr(32, p0=1/2)

## ------------------------------------------------------------------------
discrim(17, 32, method = "duotrio", conf.level = 0.9)

## ------------------------------------------------------------------------
discrim(58, 95, method = "duotrio", conf.level = 0.9)

## ------------------------------------------------------------------------
(HA <- 40/50)
(FA <- 20/50)
qnorm(HA) - qnorm(FA)

## ------------------------------------------------------------------------
pc_unb <- function(HA, FA){ pnorm((qnorm(HA) - qnorm(FA))/2)}
pcABX_IO <- function(d){pnorm(d/sqrt(2))*pnorm(d/2) + pnorm(-d/sqrt(2))*pnorm(-d/2)}
pcABX_diff <- function(d){pnorm(d/sqrt(2))*pnorm(d/sqrt(6)) + pnorm(-d/sqrt(2))*pnorm(-d/sqrt(6))}

# Check:
pc_unb(4/5, 2/5)
pcABX_IO(1.57)
pcABX_diff(1.765)

## ------------------------------------------------------------------------
inverse <-  function (f, lower = -100, upper = 100) {
  function (y) uniroot((function (x) f(x) - y), lower = lower, upper = upper)[1]
}

pcABX_IO_inverse <-  inverse(function (x) pcABX_IO(x), 0.01, 100)
pcABX_diff_inverse <-  inverse(function (x) pcABX_diff(x), 0.01, 100)

## ------------------------------------------------------------------------
pcABX_diff_inverse(pc_unb(HA, FA))
pcABX_diff_inverse(0.708)

## ------------------------------------------------------------------------
pcABX_IO_inverse(pc_unb(HA, FA))
pcABX_IO_inverse(0.708)

## ------------------------------------------------------------------------
## Doing the two-sample probit regression
y <- c(rep(0, 40), rep(1, 10), rep(0, 20), rep(1, 30))
trt <- c(rep("same", 50), rep("diff", 50))
myprobit <- glm((1-y) ~ trt, family=binomial(link="probit"))
## model summary
summary(myprobit)

## ------------------------------------------------------------------------
confint(myprobit)[2,]

## ------------------------------------------------------------------------
pcABX_diff_inverse(pnorm(0.572/2)) ## Lower CI limit
pcABX_diff_inverse(pnorm(1.632/2)) ## Upper CI limit

## ------------------------------------------------------------------------
pcABX_IO_inverse(pnorm(0.572/2)) ## Lower CI limit
pcABX_IO_inverse(pnorm(1.632/2)) ## Upper CI limit