################################################### ### Appendix A4. Dynamic factor analysis code ### ### Author: Nicole L. Michel (Nicole.L.Michel1@gmail.com) ### Modified from Case Study 4 in: Holmes, E.E., E.J. Ward, and ### M.D. Scheuerell. 2012. Analysis of multivariate time-series ### using the MARSS package. NOAA Fisheries, Northwest Fisheries ### Science Center, 2725 Montlake Blvd E., Seattle, WA 98112. ### https://cran.r-project.org/web/packages/MARSS/vignettes/UserGuide.pdf ### ### From: Michel, N.L., A.C. Smith, R.G. Clark, C.A. Morrissey, and ### K.A. Hobson. 2015. Differences in spatial synchrony and interspecific ### concordance inform guild-level population trends for aerial ### insectivorous birds. Ecography. DOI: 10.1111/ecog.01798. ################################################### ################################################### ### Load necessary packages ################################################### require(MARSS) library(xtable) ################################################### ### Load data ### Data consists of a column of years followed ### by columns with strata-level annual indices ### for all strata in the analysis region ################################################### BARSMatrixRaw <- read.csv("BARSMatrix_NW.csv") BARSMatrixUS <- BARSMatrixRaw[,2:ncol(BARSMatrixRaw)] ################################################### ### Transpose data so that years go across columns ################################################### BARSMatrixUST = t(BARSMatrixUS) N.ts = dim(BARSMatrixUST)[1] TT = dim(BARSMatrixUST)[2] ################################################### ### Select data transformation based on the proportion ### of standardized time series that are normally distributed. ################################################### transformlist <- c("sqrt","cubrt","log","none") pnormvec <- vector(mode="numeric", length=4) i <- 1 for (t in transformlist){ # apply the transformation if (t=="sqrt"){ BARSMatrixUSTt <- sqrt(BARSMatrixUST) } else if (t=="cubrt") { BARSMatrixUSTt <- (BARSMatrixUST)^(1/3) } else if (t=="log") { BARSMatrixUSTt <- log(BARSMatrixUST) } else { BARSMatrixUSTt <- BARSMatrixUST } # standardize the time series (subtract mean, divide by sd) Sigma = sqrt(apply(BARSMatrixUSTt, 1, var, na.rm=TRUE)) y.bar = apply(BARSMatrixUSTt, 1, mean, na.rm=TRUE) BARSDataT = (BARSMatrixUSTt - y.bar) * (1/Sigma) rownames(BARSDataT) = rownames(BARSMatrixUSTt) # numnorm <- 0 for (r in 1:nrow(BARSDataT)) { shapp <- shapiro.test(BARSDataT[r,])$p if (shapp > 0.05) { numnorm <- numnorm+1 } } pnormvec[i] <- numnorm/nrow(BARSDataT) i <- i+1 } ################################################### ### Transform (if needed) and standardize the data ################################################### t <- which(pnormvec==max(pnormvec)) if (t==1){ BARSMatrixUSTt <- sqrt(BARSMatrixUST) } else if (t==2) { BARSMatrixUSTt <- (BARSMatrixUST)^(1/3) } else if (t==3) { BARSMatrixUSTt <- log(BARSMatrixUST) } else { BARSMatrixUSTt <- BARSMatrixUST } # standardize the time series (subtract mean, divide by sd) Sigma = sqrt(apply(BARSMatrixUSTt, 1, var, na.rm=TRUE)) y.bar = apply(BARSMatrixUSTt, 1, mean, na.rm=TRUE) BARSDataT = (BARSMatrixUSTt - y.bar) * (1/Sigma) rownames(BARSDataT) = rownames(BARSMatrixUSTt) ################################################### ### Run Dynamic Factor Analysis. ### Run models with 3 different covariance structures ### (diagonal and unequal, diagonal and equal, and ### unconstrained). For each covariance structure type, ### produce models with an increasing number of common ### trends, stopping when the AICc increases. ################################################### # Define control settings fit.cntl.list = list(minit=200, maxit=10000, MCInit=TRUE, allow.degen=FALSE, trace=1, safe=TRUE, conv.test.slope.tol=0.1) # Define variables breaking loop when AICc increases between steps breakFlag <- FALSE AICcOld <- 999999 # create data frame to store model fit statistics model.data <- data.frame() ### Diagonal and unequal covariance structure for(m in 1:15) { dfa.model = list(A="zero", R="diagonal and unequal", m=m) kemz = MARSS(BARSDataT, model=dfa.model, control=fit.cntl.list, form="dfa") save(kemz, file=paste("BARS_NW_", m, "DU.rda", sep="")) # save model to file model.data = rbind(model.data, data.frame(R="diagonal and unequal", m=m, logLik=kemz$logLik, K=kemz$num.params, AICc=kemz$AICc, Convergence = kemz$convergence, Errors = length(kemz$errors), stringsAsFactors=FALSE)) if (kemz$AICc > AICcOld) { breakFlag <- TRUE } AICcOld <- kemz$AICc if (breakFlag) break } # Initiate break variables breakFlag <- FALSE # define variable telling R to break loop if AICc is increasing AICcOld <- 999999 # define variable to store AICc value from previous step ### Diagonal and equal covariance structure for(m in 1:15) { dfa.model = list(A="zero", R="diagonal and equal", m=m) kemz = MARSS(BARSDataT, model=dfa.model, control=fit.cntl.list, form="dfa") save(kemz, file=paste("BARS_NW_", m, "DE.rda", sep="")) model.data = rbind(model.data, data.frame(R="diagonal and equal", m=m, logLik=kemz$logLik, K=kemz$num.params, AICc=kemz$AICc, Convergence = kemz$convergence, Errors = length(kemz$errors), stringsAsFactors=FALSE)) if (kemz$AICc > AICcOld) { breakFlag <- TRUE } AICcOld <- kemz$AICc if (breakFlag) break } # Initiate break variables breakFlag <- FALSE # define variable telling R to break loop if AICc is increasing AICcOld <- 999999 # define variable to store AICc value from previous step ### Unconstrained covariance structure for(m in 1:15) { dfa.model = list(A="zero", R="unconstrained", m=m) kemz = MARSS(BARSDataT, model=dfa.model, control=fit.cntl.list, form="dfa") save(kemz, file=paste("BARS_NW_", m, "UN.rda", sep="")) model.data = rbind(model.data, data.frame(R="unconstrained", m=m, logLik=kemz$logLik, K=kemz$num.params, AICc=kemz$AICc, Convergence = kemz$convergence, Errors = length(kemz$errors), stringsAsFactors=FALSE)) if (kemz$AICc > AICcOld) { breakFlag <- TRUE } AICcOld <- kemz$AICc if (breakFlag) break } ################################################### ### Create and output model summary table ################################################### model.tbl.tot = model.data[order(model.data$AICc),-11] model.tbl.tot$delta.AICc = model.tbl.tot$AICc - min(model.tbl.tot$AICc) wt = exp(-0.5*model.tbl.tot$delta.AICc) model.tbl.tot$Ak.wt = wt/sum(wt) model.tbl.tot$Ak.wt.cum = cumsum(model.tbl.tot$Ak.wt) write.csv(model.tbl.tot, file="BARS_NW_ModelTable.csv")