# ============================================================
# ============================================================
# Parameter definitions:
# gencsv: gene csv file with genotype data
# phecsv: phenotype csv file with trait measurements
# kkcsv: kinship matrix csv file for mixed models (Mode A)
# zcsv : random-effect design csv (Mode B, n x q)  [NEW]
# covcsv: fixed covariate csv (optional)           [NEW]
#
# phename: name of the trait column to be analyzed
#   - For binary trait type: use column with 0/1 values (e.g., "Binary")
#   - For binomial trait type: use success count column (e.g., "Event")
#   - For poisson trait type: use count data column (e.g., "Poisson")
#   - For ordinal trait type: use column with ordered categories (e.g., "Ordinal")
# phetype: trait distribution type: "binary", "poisson", or "binomial"
# phenameother: secondary trait column (required only for binomial trait type)
#   - For binomial trait type: specify the trial count column (e.g., "Trial")
#   - Not needed if phetype is "binary" or "poisson"
# outpath: output directory for result files
# h2: logical (TRUE/FALSE) - whether to calculate heritability and other genetic parameters for each locus
#
# [NEW for PSR / PSR_ORDINAL (modified versions)]
# use_fixed   : TRUE/FALSE, whether to include covariates from covcsv
# use_kinship : TRUE/FALSE, Mode A (random covariance from K)
# use_randomZ : TRUE/FALSE, Mode B (random covariance from Z Z')
#
# IMPORTANT:
#   - Your modified functions do NO alignment. Row order in phe/cov/Z must match genotype individual order.
#   - Mode selection must be exclusive:
#       Mode A: use_kinship=TRUE,  use_randomZ=FALSE, kkcsv provided
#       Mode B: use_kinship=FALSE, use_randomZ=TRUE,  zcsv  provided
#
# ============================================================
# Minimal format examples (cov.csv / Z.csv)
# ============================================================
# cov.csv (minimal):
#   - at least 2 columns
#   - ID column is recommended (name "ID"), will be dropped automatically
#   - numeric/integer allowed; character will be converted to factor and dummy-coded
# Example header:
#   ID,location,sex,batch,OTHER
#
# Z.csv (minimal):
#   - n rows (same as individuals), q columns (random features)
#   - optional ID column named "ID" (will be dropped automatically)
# Example header:
#   ID,z1,z2,z3
#
# ============================================================


# ------------------------------------------------------------
# Generalized linear models for "binary", "binomial", and "poisson" trait
# (UNCHANGED - keep if you have GLM() implemented)
# ------------------------------------------------------------
GLM(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  phename = "Event",
  phetype = "binomial",  # or "binary" or "poisson"
  phenameother = "Trial", # If phetype is "binary" or "poisson", then the phenameother parameter is not needed
  outpath = "GLM_Binomial"
)

# Generalized linear models for "Ordinal" trait
GLM_ORDINAL(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  phename = "Ordinal",
  outpath = "GLM_Ordinal"
)

# ------------------------------------------------------------
# Generalized linear mixed model for "binary", "binomial", and "poisson" trait
# (UNCHANGED - keep if you have GLMM() implemented)
# ------------------------------------------------------------
GLMM(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  phename = "Event",
  phetype = "binomial",
  phenameother = "Trial",
  outpath = "GLMM_Binomial",
  h2 = FALSE  # or "TRUE"
)

# Generalized linear mixed model for "Ordinal" trait
GLMM_ORDINAL(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  phename = "Ordinal",
  outpath = "GLMM_Ordinal",
  h2 = FALSE
)

# ------------------------------------------------------------
# Population parameters previously determined for "binary", "binomial", and "poisson" trait
# ------------------------------------------------------------
P3D(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  phename = "Event",
  phetype = "binomial",
  phenameother = "Trial",
  outpath = "result_Binomial",
  h2 = FALSE
)

# Population parameters previously determined for "Ordinal" trait
P3D_ORDINAL(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  phename = "Ordinal",
  outpath = "P3D_Ordinal",
  h2 = FALSE
)

# ------------------------------------------------------------
# The pseudo response method (PSR) for "binary", "binomial", and "poisson" trait
# [CHANGED] Original PSR(...) -> Your modified PSR(...) (supports Mode A/B + cov)
# ------------------------------------------------------------

# ===== Mode A (Kinship K), no covariates =====
PSR(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  zcsv   = NULL,
  covcsv = NULL,

  phename = "Event",
  phetype = "binomial",
  phenameother = "Trial",

  outpath = "PSR_A1_K_noCov_Binomial",
  h2 = FALSE,

  use_fixed   = FALSE,
  use_kinship = TRUE,
  use_randomZ = FALSE
)

# ===== Mode A (Kinship K), with covariates =====
PSR(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  zcsv   = NULL,
  covcsv = "data/cov.csv",

  phename = "Event",
  phetype = "binomial",
  phenameother = "Trial",

  outpath = "PSR_A2_K_withCov_Binomial",
  h2 = FALSE,

  use_fixed   = TRUE,
  use_kinship = TRUE,
  use_randomZ = FALSE
)

# ===== Mode B (Random Z), no covariates =====
PSR(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = NULL,
  zcsv   = "data/Z.csv",
  covcsv = NULL,

  phename = "Event",
  phetype = "binomial",
  phenameother = "Trial",

  outpath = "PSR_B1_Z_noCov_Binomial",
  h2 = FALSE,

  use_fixed   = FALSE,
  use_kinship = FALSE,
  use_randomZ = TRUE
)

# ===== Mode B (Random Z), with covariates =====
PSR(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = NULL,
  zcsv   = "data/Z.csv",
  covcsv = "data/cov.csv",

  phename = "Event",
  phetype = "binomial",
  phenameother = "Trial",

  outpath = "PSR_B2_Z_withCov_Binomial",
  h2 = TRUE,   # example: turn on locus-level H2

  use_fixed   = TRUE,
  use_kinship = FALSE,
  use_randomZ = TRUE
)

# (Optional) also show Binary / Poisson quick calls under the same PSR interface:
PSR(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  zcsv   = NULL,
  covcsv = NULL,

  phename = "Binary",
  phetype = "binary",
  phenameother = NULL,

  outpath = "PSR_A1_K_noCov_Binary",
  h2 = FALSE,

  use_fixed   = FALSE,
  use_kinship = TRUE,
  use_randomZ = FALSE
)

PSR(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  zcsv   = NULL,
  covcsv = NULL,

  phename = "Poisson",
  phetype = "poisson",
  phenameother = NULL,

  outpath = "PSR_A1_K_noCov_Poisson",
  h2 = FALSE,

  use_fixed   = FALSE,
  use_kinship = TRUE,
  use_randomZ = FALSE
)

# ------------------------------------------------------------
# The pseudo response method (PSR) for "Ordinal" trait
# [CHANGED] 
#   PSR_ORDINAL(gencsv,phecsv, kkcsv=NULL, zcsv=NULL, covcsv=NULL,
#               phename,outpath,h2, use_fixed, use_kinship, use_randomZ)
# ------------------------------------------------------------

# ===== Ordinal, Mode A1 (Kinship K), no covariates =====
PSR_ORDINAL(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  zcsv   = NULL,
  covcsv = NULL,

  phename = "Ordinal",
  outpath = "PSR_Ordinal_A1_K_noCov",
  h2 = FALSE,

  use_fixed   = FALSE,
  use_kinship = TRUE,
  use_randomZ = FALSE
)

# ===== Ordinal, Mode A2 (Kinship K), with covariates =====
PSR_ORDINAL(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = "data/kk.csv",
  zcsv   = NULL,
  covcsv = "data/cov.csv",

  phename = "Ordinal",
  outpath = "PSR_Ordinal_A2_K_withCov",
  h2 = FALSE,

  use_fixed   = TRUE,
  use_kinship = TRUE,
  use_randomZ = FALSE
)

# ===== Ordinal, Mode B (Random Z), no covariates =====
PSR_ORDINAL(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = NULL,
  zcsv   = "data/Z.csv",
  covcsv = NULL,

  phename = "Ordinal",
  outpath = "PSR_Ordinal_B1_Z_noCov",
  h2 = TRUE,

  use_fixed   = FALSE,
  use_kinship = FALSE,
  use_randomZ = TRUE
)

# ===== Ordinal, Mode B (Random Z), with covariates =====
PSR_ORDINAL(
  gencsv = "data/IMF2-Genotypes-N.csv",
  phecsv = "data/IMF2-Phenotypes.csv",
  kkcsv  = NULL,
  zcsv   = "data/Z.csv",
  covcsv = "data/cov.csv",

  phename = "Ordinal",
  outpath = "PSR_Ordinal_B2_Z_withCov",
  h2 = TRUE,

  use_fixed   = TRUE,
  use_kinship = FALSE,
  use_randomZ = TRUE
)