tikz in LaTeX and Structural Equation Modeling

Standard

During grad school, I attended an ESA Workshop on Structural Equation Modeling (SEM) let by Jim Grace. The approach allows for multivariate analysis with multiple predictors, multiple response variables, and latent variables. Up until now, my research never required using the method and I never bought the software he recommended at the time because the GUI program recommended by Grace was too expensive for my limited needs.

Recently, I had a need to use SEM at work. We had two response variables: environmental DNA (eDNA) and the ash-free dry weight of an aquatic organism (AFDW). Both were predicted by multiple environmental variables and AFDW predicted eDNA. A perfect problem for SEM.

To refresh myself of SEM, I revisited Grace’s work. I discovered that he maintains an excellent tutorial about SEM. The pages provide a nice introduction, as does his (slightly outdated) book, his classic book, and a recent Ecoshephere article.

However, I did not have a nice way to plot my results. I did not want to use a WYSIWYG tool like Inkscape or Power Point. But I remembered the tikz package in LaTeX. Here’s the figure I created:

Example SEM plot.

I created the figure using this LaTeX code:

 \documentclass{article}

 \usepackage[paperheight =11.3cm, paperwidth =9.5cm, margin = 0.1cm]{geometry} \usepackage{tikz} \usetikzlibrary{arrows} \usetikzlibrary{positioning} \begin{document} \pagenumbering{gobble} \begin{tikzpicture}[ -> , >=stealth',auto,node distance=3.5cm, thick,main node/.style={rectangle,draw, font=\sffamily}] \node[main node] (1) {Lake}; \node[main node] (2) [below of=1] {Depth}; \node[main node] (3) [below of=2] {Non-habitat}; \node[main node] (4) [below of=3] {Habitat}; \node[main node] (6) [below right of=2, align = center] {AFDW\\ $$r^2 = 0.223$$}; \node[main node] (7) [right of=6, align = center] {eDNA\\ $$r^2 = 0.384$$}; \path[every node/.style={font=\sffamily\small}] (1) edge node [above = 40pt] {\textbf{0.497}} (6) (2) edge node [left = 10pt] {\textbf{-0.370}} (6) (3) edge node [above] {0.094} (6) (4) edge node [left = 10pt] {0.116} (6) (1) edge[bend left] node [above = 10 pt] {\textbf{0.385}} (7) (2) edge[bend left] node [above = 5pt ] {0.197} (7) (3) edge[bend right] node [above = 0pt] {-0.298} (7) (4) edge[bend right] node [below = 5pt] {0.204} (7) (6) edge node [ ] {-0.180} (7); \end{tikzpicture} 

\end{document} 

Using RcppArmadillo for a matrix population model

Standard

I’ve had a busy winter and spring with training for the Birkie and then recovering from the Birkie, hence the few posts. One of the things I’ve been doing is teaching myself the Rcpp package in R. This package lets people easily use C++ code within R. In this post, I demonstrate how I used Rcpp and specifically the RcppArmadillo package to create a population model.

Matrix models are popular in ecology. These models are series of difference equations (i.e., discrete time). I was interested in coding a simple example for a two life-stage species with the projection matrix R.

To code this up in R, I use the following code

popModelR A = matrix(c(0, 2, 0.25, 0.5), nrow = 2, byrow = TRUE),
P0 = c(0, 50)){
P = matrix(0, nrow = dim(A)[1], ncol = nYears + 1)
P[, 1] = P0
for( t in 1:nYears){
P[ , t + 1] = A  %*% P[, t]
}
return(P)
}

nYears = 10
A = matrix(c(0, 2, 0.25, 0.5),
nrow = 2, byrow = TRUE)
P0 = matrix(c(0, 50), nrow = 2)
P0
popModelROut <- popModelR(nYears = 10,
A = A,
P0 = P0)

Obviously, this simple model runs fast, but how would one code this with Rcpp? In order to get matrix operators, I needed to use the RcppArmadillo package, so my code looks like this:

library(“inline”)
library(“Rcpp”)

src1 <- ‘
int nYearsX = Rcpp::as<int>(nYears);
arma::mat P0X = Rcpp::as<arma::mat>(P0);
arma::mat AX  = Rcpp::as<arma::mat>(A);
arma::mat PX(AX.n_cols, nYearsX + 1);
PX.col(0) = P0X;

for(int t = 0; t < nYearsX; t++) {
PX.col(t + 1) =  AX * PX.col(t);
}

return Rcpp::wrap(PX);

popModelRcpp <- cxxfunction(signature(nYears = “integer”,
A = “matrix”,
P0 = “matrix”),

popModelRcpp(nYears, A, P0)

Now, to compare the two functions, I use the benchmark package and run the model for 100 simulated years:

library(rbenchmark)
nYears = 100
res popModelR(nYears, A, P0),
columns = c(“test”, “replications”, “elapsed”,
“relative”, “user.self”, “sys.self”),
order = “relative”,
replications = 1000)
print(res)

test replications elapsed relative user.self sys.self
1 popModelRcpp(nYears, A, P0) 1000 0.02 1.00 0.00 0.00
2 popModelR(nYears, A, P0) 1000 0.53 29.61 0.64 0.00

The Rcpp code is almost 30 times quicker than the base code in R!