23  Working with Nodal Attributes in Network Models

This tutorial will show how to estimate a dynamic network model for a two-group network: instead of one large group of homogeneous nodes, now we’ll have two groups that can differ based on their behavior or biology. Here, we will motivate this with a heterosexual network in which the groups may represent females and males in the population.

To get started, load the EpiModel library.

library(EpiModel)

We also recommend that you clear your environment to avoid any leftover objects from the last tutorials or labs.

rm(list = ls())

For this particular tutorial, we will set the seed to ensure reproducible results.

set.seed(12345)

23.1 Network Model

For the network model parameterization, we will review how to fit a two-group model with a differential degree distribution by group.

23.1.1 Initialization

The number in each group will be 250 nodes.

num.g1 <- num.g2 <- 250
nw <- network_initialize(n = num.g1 + num.g2)
nw
 Network attributes:
  vertices = 500 
  directed = FALSE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 0 
    missing edges= 0 
    non-missing edges= 0 

 Vertex attribute names: 
    vertex.names 

No edge attributes

We use the following function to define a nodal attribute, group, and then set it on the network. We can print out the network again to see that it has been added to the list of vertex attributes. To use group as a special attribute (see paragraph below), the values of group must be 1 and 2 only.

group <- rep(1:2, times = c(num.g1, num.g2))
nw <- set_vertex_attribute(nw, attrname = "group", value = group)
nw
 Network attributes:
  vertices = 500 
  directed = FALSE 
  hyper = FALSE 
  loops = FALSE 
  multiple = FALSE 
  bipartite = FALSE 
  total edges= 0 
    missing edges= 0 
    non-missing edges= 0 

 Vertex attribute names: 
    group vertex.names 

No edge attributes

We can use get_vertex_attribute to extract the same attribute from the network object:

get_vertex_attribute(nw, "group")
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2
[260] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[297] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[334] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[371] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[408] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[445] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[482] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

In core EpiModel, the group attribute has a special role for these built-in models that we are exploring in this course. We can, and will, use other nodal attributes on the network to model elements like age or race within the network structure. But group is special because it allows for easy parameterization of heterogeneous epidemic parameters (e.g., group-specific recovery rates). Epidemic parameters stratified by other attributes can also be arbitrarily added, but this requires some additional programming of the epidemic modules in EpiModel. We’ll get to that at the end of Day 4 and spend most of Day 5 on this.

23.1.2 Degree Distributions

Next, the group-specific degree distributions are specified as the fractions of each group that have 0, 1, 2, and 3 or more edges at any one time (momentary degree). Within this two-group framework, the first group will represent women and the second men. In our hypothetical empirical data, women exhibit less concurrency than men, but also fewer isolates. The following are those fractional distributions.

deg.dist.g1 <- c(0.40, 0.55, 0.04, 0.01)
deg.dist.g2 <- c(0.54, 0.31, 0.10, 0.05)

In our model, we will use a mean degree of 0.66. How do our degree distributions compare to the distribution expected under a Poisson probability mass function with a rate equal to this the mean degree? The dpois function in R returns the probability mass for degree 0 through 2 and ppois sums the cumulative mass for degree 3+.

pois.dists <- c(dpois(0:2, lambda = 0.66), 
                ppois(2, lambda = 0.66, lower.tail = FALSE))

This barplot compares the two observed and the one estimated fractional degree distribution, each adding to 100%. The degree distribution for men more closely matches that of the Poisson distribution, but still has slightly less concurrency than expected.

par(mar = c(3,3,2,1), mgp = c(2,1,0), mfrow = c(1,1))
barplot(cbind(deg.dist.g1, deg.dist.g2, pois.dists), 
        beside = TRUE, ylim = c(0, 0.6), col = 1:4)
legend("topright", legend = paste0("deg", 0:3), 
       pch = 15, col = 1:4, 
       cex = 0.9, bg = "white")

In this model, we will represent a multi-group network in which there is purely disortative mixing (by sex) to represent heterosexual contact. With this parameterization, and either a differential group size or degree distribution, the overall number of edges expected from each group must match. The EpiModel function check_degdist_bal takes as input two group sizes and two fractional degree distributions and checks whether the total number of implied edges is the same. Note that this constraint of matching edges is relaxed when mixing is not purely disortative (because some edges occur within group); nonetheless, TERGM parameterization still does require a single edges statistic.

check_degdist_bal(num.g1, num.g2, 
                  deg.dist.g1, deg.dist.g2)
Degree Distribution Check
=============================================
        g1.dist   g1.cnt   g2.dist   g2.cnt
Deg0       0.40    100.0      0.54    135.0
Deg1       0.55    137.5      0.31     77.5
Deg2       0.04     10.0      0.10     25.0
Deg3       0.01      2.5      0.05     12.5
Edges      1.00    165.0      1.00    165.0
=============================================
** Edges balanced ** 

23.1.3 Formation and Dissolution

Since we have not specified any nodal attributes for the network other than group, the formation model will only consider aspects of the group-specific degree distribution and mixing on group. ERGM model term degree allows for stratification by a nodal attribute. The 0:1 notation within a single degree term below is a short-hand way of writing multiple degree terms (e.g., degree(0) + degree(1)). The nodematch term is a may to represent the edges on a diagonal of a mixing matrix for the specified attribute (i.e., within group mixing).

formation <- ~edges + degree(0:1, by = "group") + nodematch("group")

For the target statistics, we specify the number of edges as a function of mean degree as 0.66 (500/2 * 0.66 = 165). Or we can input the numbers directly from the check_degdist_bal output, where the Edges indicates the edges and the g1.cnt the numbers of group 1 nodes with degree 0 to 3, and g2.cnt the same for group 2 nodes. There are four target statistics needed for the degree terms: 2 for each degree term times 2 for each group. They are input by group and then by degree (so degree 0 for group 1, then degree 1 for group 1, then degree 0 for group 2, and degree 1 for group 2). In this particular model, we will not overfit by including all possible degree terms (i.e., a count of persons with degrees 2 and 3). Finally, the target statistic for nodematch with a purely disortative mixing model is 0: there are no edges between persons matched on the value of the group variable.

target.stats <- c(165, 100, 137.5, 135, 77.5, 0)

To recap the formulation of the target statistics:

  • edges is a count of edges
  • degree is a count of nodes (whith a particular degree)
  • nodematch is a count of edges (matched on a nodal attribute)

The dissolution model is parameterized the same as the prior example but with a shorter partnership duration than the first tutorial.

coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 10)
coef.diss
Dissolution Coefficients
=======================
Dissolution Model: ~offset(edges)
Target Statistics: 10
Crude Coefficient: 2.197225
Mortality/Exit Rate: 0
Adjusted Coefficient: 2.197225

23.1.4 Estimation

The netest function is used for network model fitting.

est <- netest(nw, formation, target.stats, coef.diss)
Warning: 'glpk' selected as the solver, but package 'Rglpk' is not available;
falling back to 'lpSolveAPI'. This should be fine unless the sample size and/or
the number of parameters is very big.

For the dynamic diagnostics, we simulate from the model fit. The edges have a small bias, due to the relatively short duration. It is important to consider both the relative and absolute magnitude of the bias here. For our purposes, this fit is good enough to move forward!

dx <- netdx(est, nsims = 10, nsteps = 500, ncores = 5,
            nwstats.formula = ~edges + degree(0:3, by = "group"))

Network Diagnostics
-----------------------
- Simulating 10 networks
- Calculating formation statistics
dx
EpiModel Network Diagnostics
=======================
Diagnostic Method: Dynamic
Simulations: 10
Time Steps per Sim: 500

Formation Diagnostics
----------------------- 
            Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges        165.0  158.975   -3.651  0.478 -12.601         0.748         9.995
deg0.group1  100.0  106.592    6.592  0.316  20.876         0.518         8.014
deg1.group1  137.5  130.463   -5.118  0.309 -22.774         0.636         7.890
deg2.group1     NA   10.656       NA  0.122      NA         0.261         3.114
deg3.group1     NA    1.986       NA  0.045      NA         0.184         1.361
deg0.group2  135.0  138.446    2.553  0.297  11.596         0.618         6.554
deg1.group2   77.5   75.847   -2.133  0.265  -6.242         0.863         6.984
deg2.group2     NA   26.261       NA  0.144      NA         0.397         4.435
deg3.group2     NA    7.525       NA  0.077      NA         0.224         2.501

Duration Diagnostics
----------------------- 
      Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges     10    9.966   -0.335  0.047  -0.716         0.129          0.78

Dissolution Diagnostics
----------------------- 
      Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges    0.1      0.1     0.42      0   1.216         0.001         0.024

This is confirmed when we plot the diagnostics.

plot(dx)

To compare simulations against out-of-model predictions, let’s look again at the expected degree 2 and degree 3 for each group. The simulation means are pretty close to those expectations.

check_degdist_bal(num.g1, num.g2, 
                  deg.dist.g1, deg.dist.g2)
Degree Distribution Check
=============================================
        g1.dist   g1.cnt   g2.dist   g2.cnt
Deg0       0.40    100.0      0.54    135.0
Deg1       0.55    137.5      0.31     77.5
Deg2       0.04     10.0      0.10     25.0
Deg3       0.01      2.5      0.05     12.5
Edges      1.00    165.0      1.00    165.0
=============================================
** Edges balanced ** 

23.2 Epidemic Model

For the disease simulation, we will simulate an SIR epidemic in a closed population.

23.2.1 Parameterization

The epidemic model parameters are below. In this tutorial, we will use biological parameters for the infection probability per act and recovery rates as equal to demonstrate the impact of group-specific degree distributions on outcomes. Note that the group-specific parameters for the infection probability govern the risk of infection to persons in that mode given contact with persons in the other group.

param <- param.net(inf.prob = 0.2, inf.prob.g2 = 0.2,
                   rec.rate = 0.02, rec.rate.g2 = 0.02)

At the outset, 10 people in each group will be infected. For an SIR epidemic simulation, it is necessary to specify the number recovered in each group, here starting at 0.

init <- init.net(i.num = 10, i.num.g2 = 10, 
                 r.num = 0, r.num.g2 = 0)

For control settings, we will simulate 5 epidemics over 500 time steps each. Here we use the multi-core functionality in our controls.

control <- control.net(type = "SIR", nsims = 5, nsteps = 500, ncores = 5)

23.2.2 Simulation

The model is simulated by inputting the fitted network model, the epidemic parameters, the initial conditions, and control settings.

sim <- netsim(est, param, init, control)

Printing the model object shows the inputs and outputs from the simulation object. In particular, take note of the Variables listing; the epidemic outputs for Group 1 are listed without a suffix whereas the outputs for Group 2 are listed with a .g2 suffix.

sim
EpiModel Simulation
=======================
Model class: netsim

Simulation Summary
-----------------------
Model type: SIR
No. simulations: 5
No. time steps: 500
No. NW groups: 2

Fixed Parameters
---------------------------
inf.prob = 0.2
rec.rate = 0.02
inf.prob.g2 = 0.2
rec.rate.g2 = 0.02
act.rate = 1
groups = 2

Model Output
-----------------------
Variables: s.num i.num r.num num s.num.g2 i.num.g2 
r.num.g2 num.g2 si.flow si.flow.g2 ir.flow ir.flow.g2
Networks: sim1 ... sim5
Transmissions: sim1 ... sim5

Formation Statistics
----------------------- 
                Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means)
edges            165.0  158.024   -4.228  0.626 -11.134         1.087
deg0.group1      100.0  107.013    7.013  0.441  15.909         0.591
deg1.group1      137.5  130.476   -5.108  0.398 -17.629         0.327
deg0.group2      135.0  138.878    2.873  0.418   9.269         0.666
deg1.group2       77.5   75.554   -2.511  0.432  -4.509         1.048
nodematch.group    0.0    0.000      NaN    NaN     NaN         0.000
                SD(Statistic)
edges                   9.649
deg0.group1             7.783
deg1.group1             7.746
deg0.group2             6.533
deg1.group2             7.476
nodematch.group         0.000


Duration Statistics
----------------------- 
      Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges     10    10.03    0.303  0.064   0.471         0.138         0.772

Dissolution Statistics
----------------------- 
      Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges    0.1      0.1   -0.005      0  -0.012         0.001         0.024

Note also that the print output includes the names of the individual modules. For these two-group models, EpiModel automatically selects which modules and functions to use. For extension modules (that do not have have a type input), you will need to specify this.

23.2.3 Analysis

Similar to the first tutorial, plotting the netsim object shows the prevalence for each disease state in the model over time.

par(mfrow = c(1, 1))
plot(sim, main = "Disease State Sizes")

This plot more clearly shows prevalence and incidence by mode. Interestingly, women (Group 1) have a slightly higher prevalence of disease during the main epidemic period, and more women than men end up in the recovered state. Men and women have the same mean degree but men have a higher prevalence of concurrency (degree 2+) than women. Following epidemiological and network theory, conditional on mean degree (which is the same for both groups in this model), concurrency increases the risk of transmission but not acquisition. So therefore, women bear the higher burden of male concurrency. Your course instructor, Steve Goodreau, has prepared an excellent tutorial and exercise set on this topic.

par(mfrow = c(1,2))
plot(sim, y = c("i.num", "i.num.g2"), popfrac = TRUE,
     qnts = FALSE, ylim = c(0, 0.4), legend = TRUE)
plot(sim, y = c("si.flow", "si.flow.g2"), 
     qnts = FALSE, ylim = c(0, 2.5), legend = TRUE)

It is possible to plot the color-coded static network at various time points during the simulation, as in the last tutorial. This helps to visualize how womens’ (circle) risk is dependent on their male (square) partners’ network connections.

par(mfrow = c(1, 1), mar = c(0, 0, 0, 0))
plot(sim, type = "network", col.status = TRUE, at = 50, 
     sims = "mean", shp.g2 = "square")

It is possible to add new variables to an existing model object with mutate_epi, which takes inspiration from the mutate functions in the tidyverse. With a time step unit of a week, we can calculate standardized incidence rates per 100 person-years at risk using the following approach. The mutate_epi function takes a netsim object and adds a new variable.

sim <- mutate_epi(sim, ir.g1 = (si.flow / s.num) * 100 * 52,
                       ir.g2 = (si.flow.g2 / s.num.g2) * 100 * 52)

Printing out the object, you can see that these two new variables are in the variable list.

sim
EpiModel Simulation
=======================
Model class: netsim

Simulation Summary
-----------------------
Model type: SIR
No. simulations: 5
No. time steps: 500
No. NW groups: 2

Fixed Parameters
---------------------------
inf.prob = 0.2
rec.rate = 0.02
inf.prob.g2 = 0.2
rec.rate.g2 = 0.02
act.rate = 1
groups = 2

Model Output
-----------------------
Variables: s.num i.num r.num num s.num.g2 i.num.g2 
r.num.g2 num.g2 si.flow si.flow.g2 ir.flow ir.flow.g2 ir.g1 
ir.g2
Networks: sim1 ... sim5
Transmissions: sim1 ... sim5

Formation Statistics
----------------------- 
                Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means)
edges            165.0  158.024   -4.228  0.626 -11.134         1.087
deg0.group1      100.0  107.013    7.013  0.441  15.909         0.591
deg1.group1      137.5  130.476   -5.108  0.398 -17.629         0.327
deg0.group2      135.0  138.878    2.873  0.418   9.269         0.666
deg1.group2       77.5   75.554   -2.511  0.432  -4.509         1.048
nodematch.group    0.0    0.000      NaN    NaN     NaN         0.000
                SD(Statistic)
edges                   9.649
deg0.group1             7.783
deg1.group1             7.746
deg0.group2             6.533
deg1.group2             7.476
nodematch.group         0.000


Duration Statistics
----------------------- 
      Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges     10    10.03    0.303  0.064   0.471         0.138         0.772

Dissolution Statistics
----------------------- 
      Target Sim Mean Pct Diff Sim SE Z Score SD(Sim Means) SD(Statistic)
edges    0.1      0.1   -0.005      0  -0.012         0.001         0.024

After we add the new variable, we can plot and analyze it like any other built-in variable.

par(mar = c(3,3,2,1))
plot(sim, y = c("ir.g1", "ir.g2"), legend = TRUE,
     main = "Standardized Incidence Rates by Group")

Incidence rates in particular are often highly variable. The plot by default will use a loess smoother so some of this variability may be lost. Therefore, it is worth inspecting the non-smoothed plots too.

plot(sim, y = c("ir.g1", "ir.g2"), legend = TRUE, 
     mean.smooth = FALSE, qnts = FALSE, mean.lwd = 1,
     main = "Standardized Incidence Rates by Group (Non-Smoothed)")

As another exercise, let’s extract the data into a data frame and then calculate some other relevant summary statistics. By default as.data.frame will extract the individual simulation data in each row, but by specifying out = "mean" we can see the time specific averages for all the variables across simulations. The variables in this data frame include both the original variables from netsim plus the two new ones we added with mutate_epi.

df <- as.data.frame(sim, out = "mean")
head(df, 10)
   time s.num i.num r.num num s.num.g2 i.num.g2 r.num.g2 num.g2 si.flow
1     1 240.0  10.0   0.0 250    240.0     10.0      0.0    250     NaN
2     2 238.6  11.4   0.0 250    238.2     11.8      0.0    250     1.4
3     3 237.0  13.0   0.0 250    237.6     12.2      0.2    250     1.6
4     4 235.8  13.6   0.6 250    236.4     13.4      0.2    250     1.2
5     5 235.2  13.8   1.0 250    236.0     13.6      0.4    250     0.6
6     6 235.0  13.6   1.4 250    235.4     14.0      0.6    250     0.2
7     7 234.2  14.2   1.6 250    234.4     14.4      1.2    250     0.8
8     8 233.2  15.2   1.6 250    234.4     14.0      1.6    250     1.0
9     9 232.6  15.6   1.8 250    234.2     14.0      1.8    250     0.6
10   10 231.0  17.0   2.0 250    234.0     14.2      1.8    250     1.6
   si.flow.g2 ir.flow ir.flow.g2     ir.g1     ir.g2
1         NaN     NaN        NaN       NaN       NaN
2         1.8     0.0        0.0 30.533385 39.421008
3         0.6     0.0        0.2 35.124235 13.183307
4         1.2     0.6        0.0 26.478341 26.366460
5         0.4     0.4        0.2 13.352571  8.813559
6         0.6     0.4        0.2  4.444444 13.276756
7         1.0     0.2        0.6 17.854406 22.261026
8         0.0     0.0        0.4 22.241786  0.000000
9         0.2     0.2        0.2 13.506998  4.463519
10        0.2     0.2        0.0 36.097817  4.444444

One interesting summary statistic for an SIR epidemic is the cumulative incidence. Here are the mean cumulative incidences in the first and second groups across simulations. It is necessary here to use na.rm = TRUE to calculate these statistics by removing the first row where the value is undefined.

sum(df$si.flow, na.rm = TRUE)
[1] 216.4
sum(df$si.flow.g2, na.rm = TRUE)
[1] 199.4

Here are two more that may be of interest: the time point of the peak standardized incidence rates, and also the value of the incidence rate at that time point.

which.max(df$ir.g1)
[1] 132
df$ir.g1[which.max(df$ir.g1)]
[1] 196.6157
plot(sim, y = "ir.g1",
     mean.smooth = FALSE, qnts = FALSE, mean.lwd = 1,
     main = "Time and Value of Peak Incidence in Females (Group 1)")
points(which.max(df$ir.g1), df$ir.g1[which.max(df$ir.g1)], cex = 3)

The transmission matrix again shows the individual transmission events that occurred at each time step, along with information about the infecting partner that may be relevant.

tm <- get_transmat(sim)
head(tm, 15)
# A tibble: 15 × 7
# Groups:   at, sus [15]
      at   sus   inf infDur transProb actRate finalProb
   <dbl> <int> <int>  <dbl>     <dbl>   <dbl>     <dbl>
 1     2   184   367     15       0.2       1       0.2
 2     2   222   310     22       0.2       1       0.2
 3     2   260   191     28       0.2       1       0.2
 4     2   286    65    115       0.2       1       0.2
 5     2   322    98     17       0.2       1       0.2
 6     2   476   157     55       0.2       1       0.2
 7     3   198   389      7       0.2       1       0.2
 8     4   201   294      8       0.2       1       0.2
 9     7   143   260      5       0.2       1       0.2
10     7   233   286      5       0.2       1       0.2
11     8   156   476      6       0.2       1       0.2
12    10    57   389     14       0.2       1       0.2
13    10   229   260      8       0.2       1       0.2
14    11   489    65    124       0.2       1       0.2
15    12    88   489      1       0.2       1       0.2

The infDur column shows the duration of infection (in time steps) of the infecting partner at the point of transmission. For example, an event with an infDur of 1 means that the infecting partner had just become infected themselves in the prior time step. This is the proportion of transmission that occured in that state.

mean(tm$infDur == 1)
[1] 0.05651106

One can also convert this into a phylogeny for easier visualization of these transmission chains. First, we convert the transmission matrix into a phylo class object that can then be further analyzed and visualized with the ape package. In this transmission matrix, there are multiple trees found because we seeded the infection with 20 people as initial conditions.

tmPhylo <- as.phylo.transmat(tm)
found multiple trees, returning a list of 17phylo objects

Here is a standard phylogram of the first four trees. Some initial seeds lead to many downstream infections whereas others do not. Is this a function of underlying contact networks of the seeded nodes, the biology or behavior of those nodes, or something else? This could all be tested!

par(mfrow = c(2, 2), mar = c(0, 0, 0, 0))
plot(tmPhylo[[1]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5)
plot(tmPhylo[[2]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5)
plot(tmPhylo[[3]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5)
plot(tmPhylo[[4]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5)

In phylo objects with multiple trees, each tree is stored with a named object. Note that these names correspond to the ID numbers in the transmission matrix itself. We can extract one tree like so:

names(tmPhylo)
 [1] "seed_367" "seed_310" "seed_191" "seed_65"  "seed_98"  "seed_157"
 [7] "seed_389" "seed_294" "seed_121" "seed_483" "seed_38"  "seed_465"
[13] "seed_255" "seed_187" "seed_321" "seed_409" "seed_124"
tmPhylo5 <- tmPhylo[[5]]

Then we can plot that particular tree as a phylogram and other visual types included in the ape package.

par(mfrow = c(2, 2), mar = c(0, 0, 0, 0))
plot(tmPhylo5, type = "phylogram",
     show.node.label = TRUE, root.edge = TRUE, cex = 0.5)
plot(tmPhylo5, type = "cladogram",
     show.node.label = TRUE, root.edge = TRUE, cex = 0.5)
plot(tmPhylo5, type = "fan",
     show.node.label = TRUE, root.edge = TRUE, cex = 0.5)
plot(tmPhylo5, type = "unrooted",
     show.node.label = TRUE, root.edge = TRUE, cex = 0.5)