Rotate factors to match Principal-Components Analysis

Rotate lower-triangle loadings matrix to order factors from largest to smallest variance.

Usage

rotate_pca(
  L_tf,
  x_sf = matrix(0, nrow = 0, ncol = ncol(L_tf)),
  order = c("none", "increasing", "decreasing")
)

Arguments

L_tf: Loadings matrix with dimension \(T \times F\).
x_sf: Spatial response with dimensions \(S \times F\).
order: Options for resolving label-switching via reflecting each factor to achieve a given order across dimension \(T\).

Value

List containing the rotated loadings L_tf, the inverse-rotated response matrix x_sf, and the rotation H

Examples

# Extract covariance for spatial factor analysis (too slow for CRAN)
# \donttest{
# Simulate settings
set.seed(101)
theta_xy = 0.4
n_x = n_y = 10
n_c = 3           # Number of species
n_f = 1           # Number of factors
rho = 0.8
resid_sd = 0.5

# Simulate GMRFs
R_s = exp(-theta_xy * abs(outer(1:n_x, 1:n_y, FUN="-")) )
R_ss = kronecker(X=R_s, Y=R_s)
delta_fs = mvtnorm::rmvnorm(n_c, sigma=R_ss )

# Simulate loadings for two factors
L_cf = matrix( rnorm(n_c^2), nrow=n_c )
L_cf[,seq(from=n_f+1, to=n_c)] = 0
L_cf = L_cf + resid_sd * diag(n_c)

# Simulate correlated densities
d_cs = L_cf %*% delta_fs

# Shape into longform data-frame and add error
Data = data.frame( expand.grid(species=1:n_c, x=1:n_x, y=1:n_y),
                   "var"="logn", "z"=exp(as.vector(d_cs)) )
Data$n = rnorm( n=nrow(Data), mean=Data$z, sd=1 )

# make mesh
mesh = fmesher::fm_mesh_2d( Data[,c('x','y')] )

# Specify factor model with two factors and additional independent variance with shared SD
sem = "
  # Loadings matrix
  f1 -> 1, l_1.1
  f1 -> 2, l_1.2
  f1 -> 3, l_1.3
  f2 -> 2, l_2.2
  f2 -> 3, l_3.2

  # Factor variance = 1
  f1 <-> f1, NA, 1
  f2 <-> f2, NA, 1

  # Shared residual variance
  1 <-> 1, sd, 1
  2 <-> 2, sd, 1
  3 <-> 3, sd, 1
"

# fit model
vars = c( "f1", "f2", 1:n_c )
out = tinyVAST(
  space_term = sem,
  data = Data,
  formula = n ~ 0 + factor(species),
  spatial_domain = mesh,
  variables = vars,
  space_columns = c("x","y"),
  variable_column = "species",
  time_column = "time",
  distribution_column = "dist",
  control = tinyVASTcontrol(
    extra_report = TRUE
  )
)

# Extract loadings from factor to species
L_cf = out$rep$Rho_cc[ match(1:n_c, vars), match(c("f1","f2"), vars), drop = FALSE]

# Extract spatial values for factors
omega_sf = out$internal$parlist$omega_sc[, match(c("f1","f2"), vars), drop = FALSE]

# Apply rotation
Rotation = rotate_pca( L_cf, x_sf = omega_sf )

# Extract rotated loadings
Lprime_cf = Rotation$L_tf

# Extract rotated factor values for each mesh location
Oprime_sf = Rotation$x_sf
# }