# 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