Dynamic Factor Analysis

The MARSS() function allows you to fit DFAs with the same form as MARSS(x, form="dfa"). This has a diagonal Q with 1 on the diagonal and a stochastic x₁ with mean 0 and variance of 5 (diagonal variance-covariance matrix). There are only 3 options allowed for R: diagonal and equal, diagonal and unequal, and unconstrained.

Example data

library(MARSS)
library(dplyr)
data(lakeWAplankton, package = "MARSS")
phytoplankton <- c("Cryptomonas", "Diatoms", "Greens", "Unicells", "Other.algae")
dat <- as.data.frame(lakeWAplanktonTrans) |>
  subset(Year >= 1980 & Year <= 1989) |>
  subset(select=phytoplankton) |>
  t() |>
  MARSS::zscore()

Fit models without covariates

mod.list <- list(R='unconstrained', m=1, tinitx=1)

Fit with MARSS with EM or optim and BFGS.

m1 <- MARSS(dat, model=mod.list, form='dfa', z.score=FALSE, silent = TRUE)
m2 <- MARSS(dat, model=mod.list, form='dfa', z.score=FALSE, silent = TRUE, method="BFGS")

Fit with TMB.

library(marssTMB)
m3 <- dfaTMB(dat, model=list(m=1, R='unconstrained'))
m4 <- MARSS(dat, model=mod.list, form='dfa', method="TMB", silent=TRUE)
m5 <- MARSS(dat, model=mod.list, method="BFGS_TMB", form='dfa', silent=TRUE)

Log likelihoods

name	logLik
MARSS-EM	-772.4017
MARSS-BFGS	-772.4011
dfaTMB-nlminb	-772.4011
MARSS_tmb-nlminb	-772.4011
MARSS_tmb-optim-BFGS	-783.0471

Compare parameter estimates

Add example with covariates

For form="dfa", pass in covariates with covariates=xyz. If using the default model form (not dfa), then pass in covariates with model$d or model$c.

Fit model

# use a simpler R
mod.list2 <- list(m=1, R='diagonal and unequal', tinitx=1)
# add a temperature covariate
temp <- as.data.frame(lakeWAplanktonTrans) |>
    subset(Year >= 1980 & Year <= 1989) |>
    subset(select=Temp)
covar <- t(temp) |> zscore()
m.cov1 <- MARSS(dat, model=mod.list2, form="dfa", covariates=covar, silent = TRUE, z.score = FALSE, method="TMB")

Add a 2nd covariate

TP <- as.data.frame(lakeWAplanktonTrans) |>
    subset(Year >= 1980 & Year <= 1989) |>
    subset(select=TP)
covar <- rbind(covar, t(TP)) |> zscore()
m.cov2 <- MARSS(dat, model=mod.list2, form="dfa", covariates=covar, silent = TRUE, z.score=FALSE, method="TMB")

Parameter estimates

Model with one covariate.

tidy(m.cov1)
#>                           term    estimate  std.error    conf.low     conf.up
#> 1                         Z.11  0.31418922 0.05433970  0.20768536  0.42069309
#> 2                         Z.21  0.23028343 0.05024996  0.13179531  0.32877154
#> 3                         Z.31  0.15133602 0.04841083  0.05645254  0.24621950
#> 4                         Z.41  0.56773162 0.06970745  0.43110752  0.70435571
#> 5                         Z.51 -0.05856336 0.04136120 -0.13962982  0.02250309
#> 6  R.(Cryptomonas,Cryptomonas)  0.75260003 0.10172886  0.55321512  0.95198493
#> 7          R.(Diatoms,Diatoms)  0.77087944 0.10199558  0.57097177  0.97078711
#> 8            R.(Greens,Greens)  0.71350994 0.09469766  0.52790594  0.89911394
#> 9        R.(Unicells,Unicells)  0.20849764 0.05388283  0.10288923  0.31410605
#> 10 R.(Other.algae,Other.algae)  0.68761333 0.08891885  0.51333560  0.86189107
#> 11               D.Cryptomonas  0.05877345 0.09805646 -0.13341369  0.25096058
#> 12                   D.Diatoms -0.27858516 0.09131446 -0.45755821 -0.09961210
#> 13                    D.Greens  0.51108061 0.08392952  0.34658177  0.67557945
#> 14                  D.Unicells  0.13121433 0.10950305 -0.08340770  0.34583636
#> 15               D.Other.algae  0.53925536 0.07722438  0.38789835  0.69061236

Model with two covariates.

tidy(m.cov2)
#>                           term    estimate  std.error    conf.low      conf.up
#> 1                         Z.11  0.34060027 0.05594385  0.23095234  0.450248194
#> 2                         Z.21  0.26339330 0.05139875  0.16265360  0.364132996
#> 3                         Z.31  0.18775386 0.04887535  0.09195995  0.283547781
#> 4                         Z.41  0.54020043 0.06833001  0.40627607  0.674124789
#> 5                         Z.51 -0.04027862 0.04035477 -0.11937251  0.038815272
#> 6  R.(Cryptomonas,Cryptomonas)  0.71286331 0.09804098  0.52070651  0.905020102
#> 7          R.(Diatoms,Diatoms)  0.72982361 0.09775673  0.53822393  0.921423282
#> 8            R.(Greens,Greens)  0.68092518 0.09112594  0.50232163  0.859528736
#> 9        R.(Unicells,Unicells)  0.23937755 0.05464675  0.13227188  0.346483223
#> 10 R.(Other.algae,Other.algae)  0.67572369 0.08731346  0.50459246  0.846854918
#> 11        D.(Cryptomonas,Temp) -0.07369866 0.12487871 -0.31845643  0.171059121
#> 12            D.(Diatoms,Temp) -0.40186329 0.11652894 -0.63025583 -0.173470761
#> 13             D.(Greens,Temp)  0.39488271 0.10539729  0.18830782  0.601457612
#> 14           D.(Unicells,Temp)  0.13345292 0.13402433 -0.12922993  0.396135775
#> 15        D.(Other.algae,Temp)  0.45802179 0.09598266  0.26989923  0.646144347
#> 16          D.(Cryptomonas,TP) -0.22323975 0.12563874 -0.46948715  0.023007649
#> 17              D.(Diatoms,TP) -0.20994256 0.11623684 -0.43776258  0.017877459
#> 18               D.(Greens,TP) -0.21065665 0.10654015 -0.41947150 -0.001841796
#> 19             D.(Unicells,TP)  0.01447561 0.13811710 -0.25622893  0.285180145
#> 20          D.(Other.algae,TP) -0.13975971 0.09523756 -0.32642189  0.046902476

Plot estimates

The estimated factors for one covariate model.

library(ggplot2)
autoplot(m.cov1, plot.type="xtT")

The estimated factors for two covariate model.

library(ggplot2)
autoplot(m.cov2, plot.type="xtT")