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

## ----setup--------------------------------------------------------------------
library(DistributionalDDA)

## ----load-data----------------------------------------------------------------
data("age")
str(age)

## ----explore-data-------------------------------------------------------------
# Number of observations
nrow(age$data)

# The age midpoints
age$ages

# First few rows of the data
head(age$data[, 1:7])

# Check that distributions sum to 1
summary(rowSums(age$distributions))

## ----plot-raw, fig.height=6---------------------------------------------------
# Plot all distributions with transparency, plus the mean
matplot(age$ages, t(age$distributions), type = "l", lty = 1,
        col = adjustcolor("gray", alpha.f = 0.4),
        xlab = "Age", ylab = "Proportion",
        main = "All Age Distributions with Mean")

# Add the mean distribution
meanDist <- colMeans(age$distributions)
lines(age$ages, meanDist, col = "black", lwd = 3)
legend("topright", c("Individual", "Mean"),
       col = c("gray", "black"), lwd = c(1, 3))

## ----plot-by-sex, fig.height=6------------------------------------------------
# Separate by sex
maleIdx <- which(age$data$sex == "male")
femaleIdx <- which(age$data$sex == "female")

# Plot all distributions colored by sex
matplot(age$ages, t(age$distributions[maleIdx, ]), type = "l",
        lty = 1, col = adjustcolor("blue", alpha.f = 0.3),
        ylim = c(0, 0.1), xlab = "Age", ylab = "Proportion",
        main = "Age Distributions by Sex")
matlines(age$ages, t(age$distributions[femaleIdx, ]),
         lty = 1, col = adjustcolor("red", alpha.f = 0.3))

# Add means for each sex
lines(age$ages, colMeans(age$distributions[maleIdx, ]), col = "blue", lwd = 3)
lines(age$ages, colMeans(age$distributions[femaleIdx, ]), col = "red", lwd = 3)
legend("topright", c("Male", "Female"), col = c("blue", "red"), lwd = 3)

## ----create-basis-------------------------------------------------------------
# Reference measure: Uniform on [0, 100]
rmAge <- ReferenceMeasureUniform(0, 100)

# Create the Bayes space
B2Age <- BayesSpace(rmAge)
B2Age

# Create a basis with 7 dimensions
# (this provides a good balance between flexibility and parsimony)
basisAge <- BayesSpaceBasis(B2Age, nbasis = 7)
basisAge

## ----convert-to-bdata---------------------------------------------------------
# Define the histogram breaks
# Note: The last bin (90-100) is 10 years wide, others are 5 years
breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50,
            55, 60, 65, 70, 75, 80, 85, 90, 100)

# Create density histograms from the proportions
ageHist <- BDataDensityHistograms(age$distributions,
                                   breaks = breaks,
                                   basis = basisAge)
class(ageHist)
dim(ageHist)

# Convert to continuous densities
ageDensities <- as.BDataDensities(ageHist)

# Convert to spline representation (needed for PCA)
ageSplines <- as.BDataSplines(ageDensities)
class(ageSplines)
dim(ageSplines)

## ----plot-bdata, fig.height=6-------------------------------------------------
# Plot all distributions as PDFs with transparency
plot(ageSplines, as = "PDF", col = adjustcolor("gray", alpha.f = 0.4),
     main = "All Age Distributions (as PDFs) with Mean")

# Add the mean distribution
meanAge <- mean(ageSplines)
plot(meanAge, as = "PDF", add = TRUE, col = "black", lwd = 3)
legend("topright", c("Individual", "Mean"), col = c("gray", "black"), lwd = c(1, 3))

## ----plot-bdata-sex, fig.height=6---------------------------------------------
# Plot by sex with transparency
maleSplines <- BDataSplines(unclass(ageSplines)[maleIdx, ], basis = basisAge)
femaleSplines <- BDataSplines(unclass(ageSplines)[femaleIdx, ], basis = basisAge)

plot(maleSplines, as = "PDF", col = adjustcolor("blue", alpha.f = 0.3),
     main = "Age Distributions by Sex (as PDFs)")
plot(femaleSplines, as = "PDF", col = adjustcolor("red", alpha.f = 0.3), add = TRUE)

# Add sex-specific means
plot(mean(maleSplines), as = "PDF", add = TRUE, col = "blue", lwd = 3)
plot(mean(femaleSplines), as = "PDF", add = TRUE, col = "red", lwd = 3)
legend("topright", c("Male", "Female"), col = c("blue", "red"), lwd = 3)

## ----run-pca------------------------------------------------------------------
pcaAge <- princomp(ageSplines)
pcaAge

## ----pca-summary--------------------------------------------------------------
summary(pcaAge)

## ----screeplot----------------------------------------------------------------
plot(pcaAge, type = "screeplot")

## ----variance-plot------------------------------------------------------------
plot(pcaAge, type = "variance")

## ----score-plot, fig.height=6-------------------------------------------------
# Basic score plot
plot(pcaAge, type = "scores")

## ----scores-by-sex, fig.height=6----------------------------------------------
# Create color vector
sexColors <- ifelse(age$data$sex == "male", "blue", "red")

# Score plot colored by sex
plot(pcaAge, type = "scores", col = sexColors)
legend("topright", c("Male", "Female"), col = c("blue", "red"), pch = 19)

## ----scores-with-labels, fig.height=6-----------------------------------------
# Use row names as labels (commune codes)
plot(pcaAge, type = "scores", labels = TRUE, col = sexColors, cex = 0.7)

## ----loadings-plot, fig.height=6----------------------------------------------
# Use colorblind-friendly colors: black (mean), blue (PC1), orange (PC2)
cbPalette <- c("black", "blue", "orange")
plot(pcaAge, type = "loadings", npcs = 2, col = cbPalette)

## ----loadings-scaled, fig.height=6--------------------------------------------
# Show +/- 2 standard deviations from mean
plot(pcaAge, type = "loadings", npcs = 2, scale.sdev = 2, col = cbPalette)

## ----loadings-clr, fig.height=6-----------------------------------------------
# Centered log-ratio (CLR) representation
plot(pcaAge, type = "loadings", npcs = 2, as = "clr", col = cbPalette)

## ----covariance, fig.height=6-------------------------------------------------
plot(pcaAge, type = "covariance")

## ----reconstruction, fig.height=6---------------------------------------------
# Reconstruct using only 2 PCs with transparency
sexColorsTransp <- ifelse(age$data$sex == "male",
                          adjustcolor("blue", alpha.f = 0.3),
                          adjustcolor("red", alpha.f = 0.3))
plot(pcaAge, type = "reconstruction", npcs = 2, col = sexColorsTransp)

## ----boxplot------------------------------------------------------------------
plot(pcaAge, type = "boxplot", npcs = 4)

## ----biplot, fig.height=6-----------------------------------------------------
biplot(pcaAge)

## ----grouped-pca--------------------------------------------------------------
pcaGrouped <- princomp(ageSplines, group = age$data$sex)
pcaGrouped

## ----grouped-screeplot, fig.width=10------------------------------------------
plot(pcaGrouped, type = "screeplot")

## ----grouped-loadings, fig.width=10-------------------------------------------
plot(pcaGrouped, type = "loadings", npcs = 2)

## ----grouped-scores, fig.width=10---------------------------------------------
plot(pcaGrouped, type = "scores")

## ----multivariate-pca---------------------------------------------------------
pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE)
pcaMulti

## ----multi-screeplot----------------------------------------------------------
plot(pcaMulti, type = "screeplot")

## ----multi-scores-------------------------------------------------------------
plot(pcaMulti, type = "scores")

## ----multi-eigenfunctions, fig.width=10, fig.height=6-------------------------
plot(pcaMulti, type = "eigenfunctions", npcs = 2)

## ----predict-example----------------------------------------------------------
# Suppose we have the first 5 observations as "training" data
# and want to project observations 6-10 as "new" data

# Create a subset as new BDataSplines (subsetting keeps the class)
newAgeSplines <- BDataSplines(unclass(ageSplines)[6:10, ], basis = basisAge)

# The pcaAge model was fit on all data, but let's demonstrate predict():
newScores <- predict(pcaAge, newAgeSplines)
print(newScores[, 1:3])  # Show first 3 PCs

# Compare with the actual scores
actualScores <- pcaAge$scores[6:10, 1:3]
print(actualScores)

# They should be identical (since we used the same data)

## ----workflow-summary, eval=FALSE---------------------------------------------
# # 1. Load data
# data("age")
# 
# # 2. Create Bayes space and basis
# rmAge <- ReferenceMeasureUniform(0, 100)
# basisAge <- BayesSpaceBasis(BayesSpace(rmAge), nbasis = 7)
# 
# # 3. Convert to BData (adapt based on your data format)
# breaks <- c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50,
#             55, 60, 65, 70, 75, 80, 85, 90, 100)
# ageSplines <-
#     BDataDensityHistograms(age$distributions,
#                            breaks = breaks,
#                            basis = basisAge) |>
#   as.BataDensities() |>
#   as.BDataSplines() # pipe |> instead of nesting
# 
# # 4. Run PCA
# pcaResult <- princomp(ageSplines)
# 
# # 5. Explore results
# summary(pcaResult)                              # Variance explained
# plot(pcaResult, type = "screeplot")             # How many components?
# plot(pcaResult, type = "loadings", npcs = 2)    # What do components mean?
# plot(pcaResult, type = "scores", col = myColors) # Where are my observations?
# 
# # 6. Optional: Separate grouped analysis (compare groups)
# pcaGrouped <- princomp(ageSplines, group = age$data$sex)
# plot(pcaGrouped, type = "loadings", npcs = 2)   # Compare groups
# 
# # 7. Optional: Multivariate SFPCA (joint analysis, requires paired data)
# pcaMulti <- princomp(ageSplines, group = age$data$sex, multivariate = TRUE)
# plot(pcaMulti, type = "eigenfunctions", npcs = 2)  # Cross-group covariation

## ----session-info-------------------------------------------------------------
sessionInfo()