The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Nonlinear dynamics

James Thorson

Nonlinearity

dsem can be specified to estimate some types of nonlinearity. Here, we demonstrate an approximation to the exponential function using Lotka-Volterra dynamics.

To show this, we predict sea surface temperature from Departure Bay based upon the Pacific Decadal Oscillation:

library(dsem)
dataset = c("hare_lynx", "paramesium_didinium" )[2]

# Load data
if( dataset == "paramesium_didinium" ){
  data(paramesium_didinium)
  orig_dat = data.frame( 
    X = paramesium_didinium[,'paramecium'] / 100,
    Y = paramesium_didinium[,'didinium'] / 100 
  )
}
if( dataset == "hare_lynx" ){
  data(hare_lynx)
  orig_dat = data.frame( X = hare_lynx$hares / 10000, Y = hare_lynx$lynx / 10000 )
}

# Format
dat = full_dat = cbind(
  logX = log(orig_dat$X), logY = log(orig_dat$Y),
  X = NA, Y = NA,
  logX1 = NA, logY1 = NA,
  logX2 = NA, logY2 = NA,
  logX3 = NA, logY3 = NA,
  ones = 1
)

# Center variables for numerical stability
mean_j = colMeans( dat[,1:2], na.rm = TRUE )
dat[,1:2] = sweep( dat[,1:2], FUN = "-", MARGIN = 2, STATS = mean_j )

We then define a MDSEM:

sem = "
  # Main interactions
  logX -> logX, 1, NA, 1
  ones -> logX, 0, alpha
  Y -> logX, 1, beta, -0.1

  # Form X \approx exp(logX)
  ones -> X, 0, NA, 1
  logX -> logX1, 0, NA, 1
  logX1 -> X, 0, NA, 1
  logX1 -> logX2, 0, logX
  logX2 -> X, 0, NA, 0.5
  logX2 -> logX3, 0, logX
  logX3 -> X, 0, NA, 0.166

  # Variances
  X <-> X, 0, NA, 0
  logX <-> logX, 0, sd_logX
  logX1 <-> logX1, 0, NA, 0
  logX2 <-> logX2, 0, NA, 0
  logX3 <-> logX3, 0, NA, 0

  # Main interactions
  logY -> logY, 1, NA, 1
  X -> logY, 1, gamma
  ones -> logY, 0, delta, -0.1

  # Form Y \approx exp(logY)
  ones -> Y, 0, NA, 1
  logY -> logY1, 0, NA, 1
  logY1 -> Y, 0, NA, 1
  logY1 -> logY2, 0, logY
  logY2 -> Y, 0, NA, 0.5
  logY2 -> logY3, 0, logY
  logY3 -> Y, 0, NA, 0.166

  # Variances
  Y <-> Y, 0, NA, 0
  logY <-> logY, 0, sd_logY
  logY1 <-> logY1, 0, NA, 0
  logY2 <-> logY2, 0, NA, 0
  logY3 <-> logY3, 0, NA, 0

  # Dummy constant
  ones <-> ones, 0, NA, 0.001
  ones -> ones, 1, NA, 1
"

We then fit this without estimating any mu parameters:

fit = dsem(
  tsdata = ts(dat),
  sem = sem,
  estimate_mu = vector(), 
  estimate_delta0 = FALSE,
  control = dsem_control(
    quiet = TRUE
  )
)

We can also visualize the estimated graph

library(igraph)
library(ggraph)

g = make_empty_graph(15)
V(g)$name = c("ones", "lnx[t+1]", "lnx", "z1", "lnx^2", "z2", "lnx^3", "x", "y", "lny^3", "z3", "lny^2", "z4", "lny", "lny[t+1]" )
V(g)$shape = c( "square", "circle", "diamond")[c(1,1,1,3,2,3,2,2,2,2,3,2,3,2,1)]
V(g)$label = c("ones", "lnx[t+1]", "lnx[t]", "", "lnx[t]^2", "", "lnx[t]^3", "tilde(x)[t]", "tilde(y)[t]", "lny[t]^3", "", "lny[t]^2", "", "lny[t]", "lny[t+1]" )

#
g <- add_edges(g, c("ones", "lnx[t+1]"), attr = list(label = "alpha", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lnx", "lnx[t+1]"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lnx", "lnx^2"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lnx^2", "lnx^3"), attr = list(label = "", type = "solid", col = "black", curve = 0))
#
g <- add_edges(g, c("ones", "x"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lnx", "x"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lnx^2", "x"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lnx^3", "x"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("x", "lny[t+1]"), attr = list(label = "gamma", type = "solid", col = "black", curve = 0))
#
g <- add_edges(g, c("lnx", "z1"), attr = list(label = "", type = "solid", col = "black", curve = 1))
g <- add_edges(g, c("lnx", "z2"), attr = list(label = "", type = "solid", col = "black", curve = 1))
#
g <- add_edges(g, c("ones", "lny[t+1]"), attr = list(label = "delta", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lny", "lny[t+1]"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lny", "lny^2"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lny^2", "lny^3"), attr = list(label = "", type = "solid", col = "black", curve = 0))
#
g <- add_edges(g, c("ones", "y"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lny", "y"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lny^2", "y"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("lny^3", "y"), attr = list(label = "", type = "solid", col = "black", curve = 0))
g <- add_edges(g, c("y", "lnx[t+1]"), attr = list(label = "beta", type = "solid", col = "black", curve = 0))
#
g <- add_edges(g, c("lny", "z4"), attr = list(label = "", type = "solid", col = "black", curve = -1))
g <- add_edges(g, c("lny", "z3"), attr = list(label = "", type = "solid", col = "black", curve = -1))

pos = function(i){
  rad = i/17*2*pi - 0.25*2*pi
  return( c(x=sin(rad),y = cos(rad)) )
}
loc_nodes = t(sapply(c(0,2:14,15), pos))
# Manual fixes
loc_nodes[4,] = loc_nodes[4,] + c(0, -0.07)
loc_nodes[6,] = loc_nodes[6,] + c(-0.05, -0.05)
loc_nodes[11,] = loc_nodes[11,] + c(-0.05, 0.05)
loc_nodes[13,] = loc_nodes[13,] + c(0, 0.07)

layout = create_layout( g, loc_nodes[,c("x","y")] )
ggraph(layout) +
  geom_edge_link2(
    arrow = arrow(length = unit(2, "mm")),
    end_cap = ggraph::circle(7, 'mm'),
    start_cap = ggraph::circle(0, 'mm'),
    aes( label = label, linetype = type, col = "black", filter = (curve==0) ),  # , edge_width = ifelse(type == "dotted", 0.8, 0.8)
    vjust = -0.2,
    hjust = 0.4,
    label_parse = TRUE
  ) +
  geom_edge_arc(
    arrow = arrow(length = unit(2, "mm")),
    end_cap = ggraph::circle(7, 'mm'),
    start_cap = ggraph::circle(0, 'mm'),
    strength = 1,
    aes( label = label, linetype = type, col = col, filter = (curve==1) ),  # , edge_width = ifelse(type == "dotted", 0.8, 0.8)
    vjust = -0.2,
    hjust = 0.4,
    label_parse = TRUE
  ) +
  geom_edge_arc(
    arrow = arrow(length = unit(2, "mm")),
    end_cap = ggraph::circle(7, 'mm'),
    start_cap = ggraph::circle(0, 'mm'),
    strength = -1,
    aes( label = label, linetype = type, col = col, filter = (curve==-1) ),  # , edge_width = ifelse(type == "dotted", 0.8, 0.8)
    vjust = -0.2,
    hjust = 0.4,
    label_parse = TRUE
  ) +
  geom_node_point(
    aes(shape = shape, size = ifelse(shape == "diamond",8,15) ),   #
    #size = 15,
    stroke = 1.2,
    color = "black",
    fill = "white"
  ) +
  geom_node_label(
    fill = "white",
    lwd = 0,
    aes(label = label),
    parse = TRUE
  ) +
  theme(panel.background = element_rect(fill = NA, color = NA)) +
  coord_cartesian( xlim = 1.2*c(-1, 1), ylim = 1.2*c(-1, 1) )  +
  scale_shape_manual(values = c("circle" = 21, "square" = 22, "diamond" = 23), guide = "none") +  # ?scale_shape
  scale_edge_linetype_manual(values = c("solid" = "solid", "dotted" = "solid"), guide = "none") +
  scale_edge_colour_manual(values = c("black" = "black", "grey" = "grey"), guide = "none") +
  scale_size( range = c(6,15), guide = "none" )

Runtime for this vignette: 1.4 secs

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.