library(QCA)
data(LR)                                         # raw Lipset data, 18 countries

# Suggested anchors: 5th, 50th and 95th percentiles (the module default)
anchors <- function(x) {
  a <- unname(quantile(x, c(0.05, 0.50, 0.95)))
  digits <- min(4, max(nchar(sub("^[^.]*\\.?", "", as.character(x)))) + 1)
  a <- round(a, digits)
  if (any(x == a[2])) a[2] <- a[2] + 10^(-digits)   # crossover never equals a case
  c(e = a[1], c = a[2], i = a[3])
}

sets <- c("DEV", "URB", "LIT", "IND", "STB", "SURV")
calibrated <- as.data.frame(lapply(sets, function(s)
  calibrate(LR[[s]], type = "fuzzy", method = "direct",
            thresholds = anchors(LR[[s]]), logistic = TRUE)),
  col.names = sets, row.names = rownames(LR))
calibrated <- round(calibrated, 2)               # membership scores, 2 decimals

anchors(LR$DEV)                                  # e = 329.4, c = 588, i = 1047

# Calibration check: cases above 0.5 and cases at exactly 0.5
data.frame(
  Set       = sets,
  In_pct    = round(100 * colMeans(calibrated > 0.5), 1),
  Ambiguous = sapply(calibrated, function(m) paste(rownames(calibrated)[m == 0.5], collapse = ", "))
)                                                # DEV: FI;  LIT: EE, IE  -> adjust c

data(LF)                                         # published fuzzy-set Lipset data
conditions <- c("DEV", "URB", "LIT", "IND", "STB")
Y <- LF$SURV

necessity_fit <- function(X, Y) {
  m <- pmin(X, Y)
  c(Consistency = sum(m) / sum(Y),
    Coverage    = sum(m) / sum(X),
    RoN         = sum(1 - X) / sum(1 - m))
}

necessity <- rbind(
  t(sapply(conditions, function(cn) necessity_fit(LF[[cn]], Y))),
  t(sapply(conditions, function(cn) necessity_fit(1 - LF[[cn]], Y))))
rownames(necessity) <- c(conditions, paste0("~", conditions))
round(necessity, 3)                              # STB: 0.920, 0.707, 0.680
                                                 # LIT: 0.991, 0.643, 0.509

# Verdict with the default thresholds 0.90 / 0.60 / 0.60
ifelse(necessity[, "Consistency"] >= 0.90 & necessity[, "Coverage"] >= 0.60 &
         necessity[, "RoN"] >= 0.60, "Necessary",
       ifelse(necessity[, "Consistency"] >= 0.90, "Possibly trivial", "Not necessary"))
                                                 # STB: "Necessary"
                                                 # LIT: "Possibly trivial" (RoN < 0.60)

pof(LF[, conditions], Y, relation = "necessity") # same consistency, coverage and RoN

TT <- truthTable(LF, outcome = "SURV", conditions = conditions,
                 incl.cut = 0.80, n.cut = 1, pri.cut = 0,
                 complete = TRUE, show.cases = TRUE, sort.by = "incl, n")
TT
# 9 of 32 configurations are observed; 23 are logical remainders.
# Row 32 (DEV*URB*LIT*IND*STB; BE, CZ, NL, UK): incl = 0.904, PRI = 0.886 -> OUT = 1
# Row 22 (DEV*~URB*LIT*~IND*STB; FI, IE):       incl = 0.804, PRI = 0.719 -> OUT = 1
# Row 24 (FR, SE):                              incl = 0.709              -> OUT = 0

# Conservative: no logical remainders
sc <- minimize(TT, include = "1", details = TRUE, show.cases = TRUE)
sc   # DEV*URB*LIT*IND*STB + DEV*~URB*LIT*~IND*STB -> SURV  (cons 0.870, cov 0.658)

# Parsimonious: every helpful remainder; $SA lists the simplifying assumptions
sp <- minimize(TT, include = "?", details = TRUE, show.cases = TRUE)
sp   # DEV*~IND + URB*STB -> SURV                          (cons 0.850, cov 0.714)
sp$SA$M1                                         # 12 remainder rows

# Intermediate: easy counterfactuals only (all conditions expected "present")
si <- minimize(TT, include = "?", dir.exp = "DEV, URB, LIT, IND, STB",
               details = TRUE, show.cases = TRUE)
si   # DEV*URB*LIT*STB + DEV*LIT*~IND*STB -> SURV          (cons 0.866, cov 0.660)
si$i.sol$C1P1$EC                                 # easy counterfactual: row 30

literals <- function(term) strsplit(term, "*", fixed = TRUE)[[1]]
intermediate <- si$i.sol$C1P1$solution[[1]]      # "DEV*URB*LIT*STB" "DEV*LIT*~IND*STB"
parsimonious <- lapply(sp$solution[[1]], literals)

# Core: literals of a parsimonious term contained in the intermediate path
lapply(setNames(intermediate, intermediate), function(path) {
  lit  <- literals(path)
  core <- unique(unlist(Filter(function(p) all(p %in% lit), parsimonious)))
  c(core = paste(core, collapse = ", "),
    peripheral = paste(setdiff(lit, core), collapse = ", "))
})
# DEV*URB*LIT*STB : core URB, STB; peripheral DEV, LIT
# DEV*LIT*~IND*STB: core DEV, ~IND; peripheral LIT, STB

X <- admisc::compute(paste(intermediate, collapse = " + "), data = LF)
Y <- LF$SURV
c(Consistency = sum(pmin(X, Y)) / sum(X),        # 0.866
  Coverage    = sum(pmin(X, Y)) / sum(Y))        # 0.660

type <- ifelse(X > 0.5 & Y > 0.5, "Typical",
        ifelse(X > 0.5 & Y < 0.5, "Deviant (consistency)",
        ifelse(X < 0.5 & Y > 0.5, "Uncovered",
        ifelse(X < 0.5 & Y < 0.5, "Irrelevant", "On the 0.5 line"))))
table(type)                                      # Typical 6, Uncovered 2, Irrelevant 10

TTn <- truthTable(LF, outcome = "~SURV", conditions = conditions,
                  incl.cut = 0.80, n.cut = 1, pri.cut = 0,
                  complete = TRUE, show.cases = TRUE, sort.by = "incl, n")
spn <- minimize(TTn, include = "?", details = TRUE, show.cases = TRUE)
spn  # ~DEV + ~STB -> ~SURV                                 (cons 0.849, cov 0.888)

# Remainders assumed to produce both SURV and ~SURV
sa_rows <- function(s) as.integer(unlist(lapply(s$SA, rownames)))
intersect(sa_rows(sp), sa_rows(spn))             # 10 12 14 16 17 21 25 29

for (cut in c(0.75, 0.80, 0.85, 0.90)) {
  tt <- truthTable(LF, outcome = "SURV", conditions = conditions,
                   incl.cut = cut, n.cut = 1, pri.cut = 0)
  s  <- minimize(tt, include = "?", dir.exp = "DEV, URB, LIT, IND, STB")
  cat(cut, sum(tt$tt$OUT == "1"), paste(s$i.sol[[1]]$solution[[1]], collapse = " + "), "\n")
}
# 0.75 and 0.80: DEV*URB*LIT*STB + DEV*LIT*~IND*STB
# 0.85 and 0.90: DEV*URB*LIT*IND*STB  (row 22 no longer passes)