## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----eval = FALSE-------------------------------------------------------------
# # install development version of package
# install.packages("devtools")
# devtools::install_github("ajbass/sffdr")

## -----------------------------------------------------------------------------
library(sffdr)
data(bmi)

# Define primary p-values and informative p-values
p <- sumstats$bmi
z <- as.matrix(sumstats[,-1])

head(sumstats)

## -----------------------------------------------------------------------------
# Create model
mpi0 <- pi0_model(z)

# Estimation of functional pi0
fpi0 <- fpi0est(p, mpi0)

## -----------------------------------------------------------------------------
# Apply sfFDR
sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0)   

# Extract results
results <- data.frame(
  p_value = p,
  fp_value = sffdr_out$fpvalues,
  fq_value = sffdr_out$fqvalues,
  flfdr = sffdr_out$flfdr
)

# Plot significance results (focusing on high-significance SNPs)
plot(sffdr_out, rng = c(0, 1e-6))

## -----------------------------------------------------------------------------
# logical vector: TRUE if SNP is LD-independent
# (In this example, all SNPs are already independent)
indep_snps <- rep(TRUE, length(p))

# Fit model using LD-independent subset
mpi0 <- pi0_model(z = z, indep_snps = indep_snps)

fpi0 <- fpi0est(p, mpi0, indep_snps = indep_snps)

# Estimate quantities for all SNPs
sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0)

## ----eval = FALSE-------------------------------------------------------------
# sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0, monotone = TRUE)

## ----eval = FALSE-------------------------------------------------------------
# # Create an artificial LD block
# LD_block <- 500:1000
# p[LD_block] <- runif(length(LD_block), min = 0.78, max = 0.82)
# 
# # Define weights (inverse LD)
# w <- rep(1.0, length(p))
# w[LD_block] <- 0.05
# 
# # Apply sfFDR with weights
# mpi0 <- pi0_model(z, weights = w)
# fpi0 <- fpi0est(p, mpi0, weights = w)
# sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0, weights = w)

## ----eval = FALSE-------------------------------------------------------------
# # Create model (can include other variables (e.g., MAF) or specify more complicated models)
# fmod <- "~ns(bfp, knots = c(0.01, 0.025, 0.05, 0.1))"
# fpi0_mod <- fpi0est(p = p,
#                     z = mpi0$zt,
#                     pi0_model = fmod)