### Simulate GWAS Data

Source: https://context7.com/cran/rrblup/llms.txt

Generates simulated genotype and phenotype data for GWAS analysis. Includes marker data, QTL effects, and phenotypes for multiple traits with environmental factors.

```R
set.seed(321)
n_lines <- 200
n_markers <- 1000
n_qtl <- 5

# Generate marker data (columns are lines for geno data frame format)
M <- matrix(0, n_markers, n_lines)
for (i in 1:n_lines) {
  M[, i] <- sample(c(-1, 0, 1), n_markers, replace = TRUE, prob = c(0.25, 0.5, 0.25))
}
colnames(M) <- paste0("Line_", 1:n_lines)

# Create genotype data frame with marker info
# Format: marker name, chromosome, position, then genotypes
geno <- data.frame(
  marker = paste0("SNP_", 1:n_markers),
  chrom = rep(1:5, each = n_markers / 5),           # 5 chromosomes
  pos = rep(1:(n_markers / 5), times = 5) * 1000,   # Position in bp
  M,
  check.names = FALSE
)

# Define QTL positions (evenly spaced across chromosomes)
qtl_positions <- c(100, 300, 500, 700, 900)  # One per chromosome
marker_effects <- rep(0, n_markers)
marker_effects[qtl_positions] <- c(1.5, -1.2, 1.8, -1.0, 1.3)  #QTL effects

# Generate phenotypes
g <- as.vector(t(M) %*% marker_effects)
h2 <- 0.5
e <- rnorm(n_lines, mean = 0, sd = sqrt((1 - h2) / h2 * var(g)))
y1 <- g + e

# Second trait with different QTL architecture
marker_effects2 <- rep(0, n_markers)
marker_effects2[c(150, 450, 750)] <- c(2.0, -1.5, 1.2)
g2 <- as.vector(t(M) %*% marker_effects2)
y2 <- g2 + rnorm(n_lines, mean = 0, sd = sqrt((1 - h2) / h2 * var(g2)))

# Create phenotype data frame
# First column must be genotype ID matching column names in geno
pheno <- data.frame(
  line = colnames(M),
  trait1 = y1,
  trait2 = y2,
  environment = factor(rep(c("Field1", "Field2"), each = n_lines / 2))
)
```

--------------------------------

### Perform Genomic Prediction with kinship.BLUP

Source: https://context7.com/cran/rrblup/llms.txt

Demonstrates generating simulated phenotypes and executing genomic prediction using RR-BLUP, Gaussian, and Exponential kernels.

```R
# True effects and phenotypes
effects <- rnorm(n_markers, sd = 0.2)
g_train <- as.vector(G_train %*% effects)
g_test <- as.vector(G_test %*% effects)
h2 <- 0.5
y <- g_train + rnorm(n_train, sd = sqrt((1 - h2) / h2 * var(g_train)))

# RR-BLUP method (equivalent to using A.mat with mixed.solve)
rr_result <- kinship.BLUP(
  y = y,
  G.train = G_train,
  G.pred = G_test,
  K.method = "RR"       # Ridge regression / additive kernel
)

cat("RR-BLUP training accuracy:", cor(g_train, rr_result$g.train), "\n")
cat("RR-BLUP prediction accuracy:", cor(g_test, rr_result$g.pred), "\n")

# Gaussian kernel method
gauss_result <- kinship.BLUP(
  y = y,
  G.train = G_train,
  G.pred = G_test,
  K.method = "GAUSS",
  n.profile = 10        # Grid points for bandwidth optimization
)

cat("\nGaussian kernel training accuracy:", cor(g_train, gauss_result$g.train), "\n")
cat("Gaussian kernel prediction accuracy:", cor(g_test, gauss_result$g.pred), "\n")
cat("Optimal kernel parameter profile:\n")
print(gauss_result$profile)

# Exponential kernel method
exp_result <- kinship.BLUP(
  y = y,
  G.train = G_train,
  G.pred = G_test,
  K.method = "EXP",
  n.profile = 10
)

cat("\nExponential kernel training accuracy:", cor(g_train, exp_result$g.train), "\n")
cat("Exponential kernel prediction accuracy:", cor(g_test, exp_result$g.pred), "\n")

# Recommended modern approach using kin.blup instead:
# A <- A.mat(rbind(G_train, G_test))
# data <- data.frame(y = y, gid = 1:n_train)
# rownames(A) <- 1:(n_train + n_test)
# result <- kin.blup(data, geno = "gid", pheno = "y", K = A)
```

--------------------------------

### mixed.solve

Source: https://context7.com/cran/rrblup/llms.txt

A general-purpose solver for mixed models of the form y = Xβ + Zu + ε, used for estimating marker effects or breeding values.

```APIDOC
## mixed.solve

### Description
Solves mixed models of the form y = Xβ + Zu + ε, where β represents fixed effects and u represents random effects with Var[u] = Kσ²u. It uses spectral decomposition to estimate variance components via ML or REML.

### Parameters
- **y** (vector) - Required - Phenotype vector.
- **Z** (matrix) - Optional - Design matrix for random effects.
- **K** (matrix) - Optional - Relationship matrix for random effects.
- **X** (matrix) - Optional - Design matrix for fixed effects.
- **SE** (boolean) - Optional - If TRUE, calculates standard errors.
- **method** (string) - Optional - Estimation method, 'ML' or 'REML'.
```

--------------------------------

### Solve Mixed Models with rrBLUP

Source: https://context7.com/cran/rrblup/llms.txt

Use `mixed.solve` for RR-BLUP (marker effects) or G-BLUP (breeding values). Supports fixed effects and standard error calculation. Requires phenotype data (y) and marker matrix (Z) or relationship matrix (K).

```r
library(rrBLUP)

# Create simulated marker data: 200 lines with 1000 SNP markers
# Markers coded as {-1, 0, 1} for {aa, Aa, AA}
set.seed(123)
n_lines <- 200
n_markers <- 1000

M <- matrix(0, n_lines, n_markers)
for (i in 1:n_lines) {
  M[i, ] <- sample(c(-1, 0, 1), n_markers, replace = TRUE, prob = c(0.25, 0.5, 0.25))
}
rownames(M) <- paste0("Line_", 1:n_lines)

# Simulate true marker effects and phenotypes
true_effects <- rnorm(n_markers, mean = 0, sd = 0.1)
genetic_values <- as.vector(M %*% true_effects)
h2 <- 0.5  # heritability
y <- genetic_values + rnorm(n_lines, mean = 0, sd = sqrt((1 - h2) / h2 * var(genetic_values)))

# Method 1: Estimate marker effects using RR-BLUP (K = I)
# Pass marker matrix as Z, K defaults to identity
marker_model <- mixed.solve(y, Z = M)

# View results
cat("Genetic variance (Vu):", marker_model$Vu, "\n")
cat("Error variance (Ve):", marker_model$Ve, "\n")
cat("Ridge parameter (lambda = Ve/Vu):", marker_model$Ve / marker_model$Vu, "\n")
cat("Log-likelihood:", marker_model$LL, "\n")
cat("Correlation with true marker effects:", cor(true_effects, marker_model$u), "\n")

# Method 2: Estimate breeding values using G-BLUP
# Use additive relationship matrix K instead of marker matrix
A <- A.mat(M)
gblup_model <- mixed.solve(y, K = A)

cat("Breeding value prediction accuracy:", cor(genetic_values, gblup_model$u), "\n")

# Method 3: Include fixed effects and get standard errors
# Create a fixed effect (e.g., environment)
env <- factor(rep(c("Env1", "Env2"), each = n_lines / 2))
X <- model.matrix(~ env)

# Solve with fixed effects and standard errors
full_model <- mixed.solve(y, Z = M, X = X, SE = TRUE, method = "REML")

cat("\nFixed effects (beta):\n")
print(full_model$beta)
cat("\nFixed effect standard errors:\n")
print(full_model$beta.SE)
cat("\nFirst 5 marker effect SEs:", head(full_model$u.SE, 5), "\n")
```

--------------------------------

### GWAS with P3D (EMMAX Equivalent)

Source: https://context7.com/cran/rrblup/llms.txt

Performs a Genome-Wide Association Study using the P3D method, which is equivalent to EMMAX. Requires pre-determined population parameters. Filters markers based on Minor Allele Frequency (MAF).

```R
cat("Running GWAS with P3D=TRUE (EMMAX)...
")
gwas_results <- GWAS(
  pheno = pheno[, c("line", "trait1")],  # Subset to single trait
  geno = geno,
  min.MAF = 0.05,       # Filter rare alleles
  P3D = TRUE,           # Use population parameters previously determined
  plot = FALSE          # Suppress plots for this example
)

# View top associations
gwas_results <- gwas_results[order(gwas_results$trait1, decreasing = TRUE), ]
cat("\nTop 10 marker associations:\n")
print(head(gwas_results, 10))

# Check if true QTL were detected
cat("\nTrue QTL positions:", qtl_positions, "\n")
cat("Were they in top 10?:", sum(qtl_positions %in% gwas_results$marker[1:10]), "of", n_qtl, "\n")
```

--------------------------------

### A.mat

Source: https://context7.com/cran/rrblup/llms.txt

Calculates the realized additive relationship matrix (A matrix) from marker data.

```APIDOC
## A.mat

### Description
Calculates the realized additive relationship matrix (A matrix) from marker data using the formula A = WW'/c. Supports missing data imputation and shrinkage estimation.

### Parameters
- **M** (matrix) - Required - Marker matrix with genotypes coded as {-1, 0, 1}.
- **min.MAF** (numeric) - Optional - Minimum minor allele frequency for marker filtering.
- **max.missing** (numeric) - Optional - Maximum proportion of missing data allowed per marker.
- **impute.method** (string) - Optional - Imputation method, 'mean' or 'EM'.
- **tol** (numeric) - Optional - Convergence tolerance for EM algorithm.
- **n.core** (integer) - Optional - Number of cores for parallel processing.
```

--------------------------------

### Estimate Shrinkage using Regression Method

Source: https://context7.com/cran/rrblup/llms.txt

Estimates the genomic relationship matrix using the regression method (Yang/Mueller), specifying the number of simulated QTL and iterations.

```R
A_shrink_reg <- A.mat(
  M[, 1:200],
  shrink = list(
    method = "REG",
    n.qtl = 50,       # Number of simulated QTL
    n.iter = 10       # Number of iterations
  ),
  n.core = 2
)
```

--------------------------------

### Calculate Additive Relationship Matrix (A.mat)

Source: https://context7.com/cran/rrblup/llms.txt

Compute the realized additive relationship matrix from marker data. Handles missing values via mean imputation or EM algorithm, and supports filtering by Minor Allele Frequency (MAF) and missing data percentage. Can utilize parallel processing.

```r
library(rrBLUP)

# Create marker data with some missing values
set.seed(456)
n_lines <- 150
n_markers <- 2000

M <- matrix(0, n_lines, n_markers)
for (i in 1:n_lines) {
  M[i, ] <- sample(c(-1, 0, 1), n_markers, replace = TRUE, prob = c(0.25, 0.5, 0.25))
}
rownames(M) <- paste0("Genotype_", 1:n_lines)
colnames(M) <- paste0("SNP_", 1:n_markers)

# Introduce 5% missing data (simulating GBS data)
missing_idx <- sample(length(M), size = round(0.05 * length(M)))
M_missing <- M
M_missing[missing_idx] <- NA

# Basic A matrix calculation (missing values imputed with mean)
A_basic <- A.mat(M_missing)

cat("Dimensions of A matrix:", dim(A_basic), "\n")
cat("Mean diagonal (should be ~1 + inbreeding):", mean(diag(A_basic)), "\n")
cat("Range of off-diagonal:", range(A_basic[upper.tri(A_basic)]), "\n")

# A matrix with filtering options
A_filtered <- A.mat(
  M_missing,
  min.MAF = 0.05,        # Remove markers with MAF < 5%
  max.missing = 0.2      # Remove markers with > 20% missing data
)

# A matrix with EM imputation for high-density GBS data
# Slower but more accurate for datasets with structured missing patterns
A_em <- A.mat(
  M_missing,
  impute.method = "EM",
  tol = 0.01,            # Convergence tolerance
  n.core = 2             # Parallel processing (Unix/Mac only)
)

```

--------------------------------

### Convert p-values and Identify Significant Markers

Source: https://context7.com/cran/rrblup/llms.txt

Converts -log10(p) values to p-values and identifies markers significant after Bonferroni correction. This is useful for downstream analysis and hypothesis testing.

```R
# Convert -log10(p) to p-values for downstream analysis
gwas_results$pvalue <- 10^(-gwas_results$trait1)
significant <- gwas_results[gwas_results$pvalue < 0.05 / n_markers, ]  # Bonferroni
cat("\nBonferroni-significant markers:", nrow(significant), "\n")
```

--------------------------------

### Gaussian Kernel BLUP for Non-Additive Effects

Source: https://context7.com/cran/rrblup/llms.txt

Applies Gaussian kernel BLUP for capturing non-additive genetic effects by using a distance matrix and specifying a grid for the scale parameter.

```R
D <- as.matrix(dist(M_all))
rownames(D) <- rownames(M_all)
colnames(D) <- rownames(M_all)

gauss_blup <- kin.blup(
  data = pheno_data,
  geno = "gid",
  pheno = "yield",
  K = D,                      # Distance matrix
  GAUSS = TRUE,               # Enable Gaussian kernel
  n.core = 2,                 # Parallel processing
  theta.seq = seq(0.1, 2, length.out = 5)  # Scale parameter grid
)

cat("\nGaussian kernel log-likelihood profile:
")
print(gauss_blup$profile)
```

--------------------------------

### Calculate Genomic Relationship Matrix with Imputation

Source: https://context7.com/cran/rrblup/llms.txt

Calculates the genomic relationship matrix (A) from marker data, with options for imputation of missing values and returning the imputed matrix.

```R
result <- A.mat(M_missing, impute.method = "mean", return.imputed = TRUE)
A_matrix <- result$A
M_imputed <- result$imputed

cat("\nImputed matrix dimensions:", dim(M_imputed), "\n")
cat("Any remaining NAs:", any(is.na(M_imputed)), "\n")
```

--------------------------------

### GWAS with Custom Kinship Matrix

Source: https://context7.com/cran/rrblup/llms.txt

Performs a GWAS using a user-defined kinship matrix (e.g., pedigree-based or marker-based). The kinship matrix is supplied via the 'K' argument.

```R
# Custom kinship matrix (e.g., pedigree-based)
K_custom <- A.mat(t(M))  # Calculate from markers
gwas_custom_k <- GWAS(
  pheno = pheno[, c("line", "trait1")],
  geno = geno,
  K = K_custom,         # Supply pre-computed kinship
  min.MAF = 0.05,
  plot = FALSE
)
```

--------------------------------

### Basic G-BLUP Prediction with kin.blup

Source: https://context7.com/cran/rrblup/llms.txt

Performs a basic Genomic Best Linear Unbiased Prediction (G-BLUP) using the kin.blup function with genotype IDs, phenotypes, and a kinship matrix.

```R
library(rrBLUP)

# Simulate a genomic selection scenario
set.seed(789)
n_training <- 150       # Phenotyped lines
n_prediction <- 50      # Unphenotyped lines for prediction
n_total <- n_training + n_prediction
n_markers <- 1000

# Generate marker data for all lines
M_all <- matrix(0, n_total, n_markers)
for (i in 1:n_total) {
  M_all[i, ] <- sample(c(-1, 0, 1), n_markers, replace = TRUE, prob = c(0.25, 0.5, 0.25))
}
rownames(M_all) <- paste0("Line_", 1:n_total)

# True genetic values for all lines
qtl_effects <- rnorm(n_markers, sd = 0.1)
g_true <- as.vector(M_all %*% qtl_effects)

# Phenotypes only for training set
h2 <- 0.6
y_train <- g_true[1:n_training] + rnorm(n_training, sd = sqrt((1 - h2) / h2 * var(g_true[1:n_training])))

# Create data frame with phenotypes
# Include environmental fixed effect and a covariate
pheno_data <- data.frame(
  gid = rownames(M_all)[1:n_training],
  yield = y_train,
  location = factor(rep(c("North", "South"), length.out = n_training)),
  plant_height = rnorm(n_training, mean = 100, sd = 10)  # covariate
)

# Calculate relationship matrix for all lines
A <- A.mat(M_all)

# Basic G-BLUP prediction
basic_gblup <- kin.blup(
  data = pheno_data,
  geno = "gid",           # Column name with genotype IDs
  pheno = "yield",        # Column name with phenotype
  K = A                   # Kinship matrix
)

cat("Genetic variance (Vg):", basic_gblup$Vg, "\n")
cat("Error variance (Ve):", basic_gblup$Ve, "\n")
cat("Heritability estimate:", basic_gblup$Vg / (basic_gblup$Vg + basic_gblup$Ve), "\n")

# Predictions for all lines (including unphenotyped)
predictions <- basic_gblup$g
cat("\nTraining accuracy:", cor(g_true[1:n_training], predictions[1:n_training]), "\n")
cat("Prediction accuracy:", cor(g_true[(n_training + 1):n_total],
                                 predictions[(n_training + 1):n_total]), "\n")
```

--------------------------------

### Legacy kinship.BLUP Function

Source: https://context7.com/cran/rrblup/llms.txt

Demonstrates the legacy `kinship.BLUP` function for genomic prediction. This function directly handles marker data and supports various kernel methods. Note: This function is superseded by `kin.blup()` and `A.mat()`.

```R
library(rrBLUP)

# Legacy function usage example
set.seed(555)
n_train <- 100
n_test <- 30
n_markers <- 500

# Generate training and test marker data
# NOTE: This function expects genotypes coded on [-1, 1] scale
G_train <- matrix(0, n_train, n_markers)
G_test <- matrix(0, n_test, n_markers)

for (i in 1:n_train) {
  G_train[i, ] <- sample(c(-1, 0, 1), n_markers, replace = TRUE, prob = c(0.25, 0.5, 0.25))
}
for (i in 1:n_test) {
  G_test[i, ] <- sample(c(-1, 0, 1), n_markers, replace = TRUE, prob = c(0.25, 0.5, 0.25))
}
```

--------------------------------

### G-BLUP with Fixed Effects, Covariates, and PEV

Source: https://context7.com/cran/rrblup/llms.txt

Performs G-BLUP including categorical fixed effects and continuous covariates, and returns prediction error variance (PEV).

```R
full_gblup <- kin.blup(
  data = pheno_data,
  geno = "gid",
  pheno = "yield",
  K = A,
  fixed = c("location"),      # Categorical fixed effect
  covariate = c("plant_height"),  # Continuous covariate
  PEV = TRUE                  # Return prediction error variance
)

cat("\nPredicted values (adjusted for fixed effects):", head(full_gblup$pred), "\n")
cat("PEV for first 5 lines:", head(full_gblup$PEV, 5), "\n")

# Calculate expected reliability from PEV
reliability <- 1 - full_gblup$PEV / (full_gblup$Vg * diag(A))
cat("Mean reliability:", mean(reliability[1:n_training], na.rm = TRUE), "\n")
```

--------------------------------

### Estimate Shrinkage using Endelman-Jannink Method

Source: https://context7.com/cran/rrblup/llms.txt

Estimates the genomic relationship matrix using the Endelman-Jannink shrinkage method, which is beneficial for prediction accuracy at low marker densities.

```R
A_shrink_ej <- A.mat(
  M[, 1:200],                      # Use subset of markers
  shrink = list(method = "EJ")     # Endelman-Jannink shrinkage
)
```

--------------------------------

### GWAS with Fixed Environmental Effect

Source: https://context7.com/cran/rrblup/llms.txt

Performs a GWAS modeling a specific environmental factor as a fixed effect. This allows for the analysis of trait variation influenced by environmental conditions.

```R
cat("\nRunning GWAS with environmental fixed effect...\n")
gwas_env <- GWAS(
  pheno = pheno[, c("line", "trait1", "environment")],
  geno = geno,
  fixed = "environment",   # Model environment as fixed effect
  min.MAF = 0.05,
  P3D = TRUE,
  plot = FALSE
)
```

--------------------------------

### Multi-trait GWAS

Source: https://context7.com/cran/rrblup/llms.txt

Conducts a Genome-Wide Association Study for multiple traits simultaneously. All non-fixed columns in the phenotype data frame are treated as traits. Supports parallel processing.

```R
cat("\nRunning multi-trait GWAS...\n")
gwas_multi <- GWAS(
  pheno = pheno[, c("line", "trait1", "trait2")],
  geno = geno,
  min.MAF = 0.05,
  P3D = TRUE,
  n.core = 2,           # Parallel processing
  plot = FALSE
)

cat("\nMulti-trait results columns:", colnames(gwas_multi), "\n")
```

--------------------------------

### GWAS with Population Structure Correction

Source: https://context7.com/cran/rrblup/llms.txt

Conducts a GWAS while correcting for population structure by including Principal Components (PCs) as fixed effects. Requires specifying the number of PCs to include.

```R
cat("\nRunning GWAS with population structure correction...\n")
gwas_pk <- GWAS(
  pheno = pheno[, c("line", "trait1")],
  geno = geno,
  n.PC = 3,             # Include first 3 PCs as fixed effects
  min.MAF = 0.05,
  P3D = TRUE,
  plot = FALSE
)
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.