The resevol package
Brad Duthie and Rose McKeon
This notebook tracks progress on the development of resevol R package, a BBSRC FAPESP Newton Funded Project to model how the spatial scale of heterogeneity in fungal isolate application and crop plant cultivation affects variation in the selection landscape for biopesticide resistance.
Project updates:
I don’t know why I didn’t think of this sooner, but there really
should be a feature that allows users to set the initial mean values of
traits. This is easy because trait means are of course independent of
the covariance structure of the traits, so there just needs to be some
specified increment to be added to each trait throughout the simulation.
I have introduced this in a new version 0.3.2.0 with the
trait_means argument, which simply takes a vector of the
means for each customised trait. By default, values are set to zero
(i.e., what they always were originally), but the argument can also take
a vector of the same length as the number of traits, with each element
position corresponding to the specific trait. That is,
trait_means[1] is the mean value of "T1",
trait_means[2] is the mean value of "T2", et
cetera. An error is given if the user attempts to specify a
trait_means vector that has a length different from the
number of traits. Just to show this in action, I will use the resevol
mine_gmatrix output mg_v1, which has the covariance
structure below.
## [,1] [,2] [,3] [,4]
## [1,] 0.9936103 -0.5106387 0.2324040 0.2367691
## [2,] -0.5106387 1.0193888 0.2573216 0.2723028
## [3,] 0.2324040 0.2573216 0.9334559 -0.4480018
## [4,] 0.2367691 0.2723028 -0.4480018 0.9689539We have four traits, which I will use to define consumption rate of
crops 1 and 2 (T1 and T2, respectively), and
the consumption rates of pesticides 1 and 2 (T3 and
T4, respectively). I will set the trait means to T1 = 1, T2
= 2, T3 = 0.2, and T4 = 0.2 and run a simulation with the code
below.
sim <- run_farm_sim(mine_output = mg_v1, N = 1000, xdim = 10, ydim = 10,
repro = "asexual", time_steps = 20, farms = 12,
crop_number = 2,
crop_init = "random",
crop_rotation_type = 2,
pesticide_number = 2,
pesticide_init = "random",
pesticide_rotation_type = 2,
food_consume = c("T1", "T2"),
pesticide_consume = c("T3", "T4"),
trait_means = c(1, 2, 0.2, 0.2),
food_needed_surv = 1,
food_needed_repr = 3,
print_last = TRUE);At the end of the simulation (after some evolution has occurred), we can look at the evolution of traits and final trait distributions to confirm that they are around what we expect. Here is the evolution of the different traits over time.
Here are the distributions of each trait, individually.
Finally, here is a pairs plot showing the trait correlations.
Overall, this is doing what we want it to. I have also written a new
test into testthat for the package. New documentation for the
trait_means argument has been written, and I will
incorporate this option into the working manuscript.
Version 0.3.0.1 has been submitted to CRAN and is passing all
checks. I am now a bit under halfway through writing a first draft
of the software note. This note be submitted as a pre-print, and for
publication, and will explain how to use the package. It will also
replace the current Get
Started article and CRAN
vignette. The structure of the publication introduces some key
background, the mine_gmatrix function and its underlying
logic, the ecology and evolution of pests in the model, and a simple
example. I think that two supporting information documents are also in
order. The first will explain the evolutionary algorithm in more detail
and basically be a version of the current
vignette with maybe some minor revisions. The second will be a more
advanced example with complex dynamics and traits, presenting
simulations like the ones shown in this journal on 26 AUG 2021. My hope is that this plus the
documentation is sufficient for anyone to start using the package.
Having submitted version 0.3.0.1, I almost immediately discovered the
need for an update that is now version 0.3.1.0. I will hold off on
submitting this one to CRAN just yet, as the update is quite minor. The
cause was that the mine_gmatrix function kept giving a
discrepancy between the stress threshold that caused it to terminate and
the actual predicted stress of the network it found. The reason for this
was that the algorithm was terminating due to finding unusually low
stress caused sampling error for the standard random normal loci. In
other words, a luck draw for random loci often caused networks to look
better at producing desired trait covariances than they actually were.
I’ve addressed this now by first using the mean stress for a whole
population of networks to test the termination criteria
(term_cri) in the mine_gmatrix function,
rather than the best stress out of all networks in the
population. But the mine_gmatrix function still returns the
network that produces the lowest found stress.
Second, I have introduced a new diagnostic function to give the distribution of stress predicted from loci.
stress_test <- function(mine_output, indivs = 1000, reps = 10){
target_VCV <- mine_output[[2]];
Loci_to_Traits <- mine_output[[5]];
loci <- mine_output[[1]][1];
Ln_stress <- rep(x = NA, times = reps);
while(reps > 0){
sim_loci <- rnorm(n = indivs * loci);
individuals <- matrix(data = sim_loci, nrow = indivs, ncol = loci);
trait_vals <- individuals %*% Loci_to_Traits;
VCV <- cov(trait_vals);
stress <- mean((target_VCV - VCV) * (target_VCV - VCV));
Ln_stress[reps] <- log(stress);
reps <- reps - 1;
}
return(Ln_stress);
}This simple function now forms part of the package itself, so users
can run the output of mine_gmatrix and get a better overall
picture on how the chosen network performs. The example used in the
documentation is below (I’ve increased the replicates using the argument
reps = 1000; the default is 10).
gmt <- matrix(data = 0, nrow = 4, ncol = 4);
diag(gmt) <- 1;
mg <- mine_gmatrix(gmatrix = gmt, loci = 4, layers = 3, indivs = 100,
npsize = 100, max_gen = 2, prnt_out = FALSE);
stresses <- stress_test(mine_output = mg, reps = 1000);The output stresses is just a vector of stresses for
some number (indivs) of initialised genomes. The histogram
below shows the distribution for the above example.
This should at least make it easier for users to decide if they do or
do not like the output of mine_gmatrix for use in a new set
of simulations. The above example shows stress that is quite consistent,
but if we let the simulation run for longer to find a better solution,
the stress variance seems to increase (note that now
max_gen = 400 below).
gmt <- matrix(data = 0, nrow = 4, ncol = 4);
diag(gmt) <- 1;
mg <- mine_gmatrix(gmatrix = gmt, loci = 4, layers = 3, indivs = 100,
npsize = 100, max_gen = 400, prnt_out = FALSE);
stresses <- stress_test(mine_output = mg, reps = 1000);
hist(stresses, xlab = "Stress values for initialised genomes",
main = "", breaks = 20);
It might be good to recommend this kind of diagnostic for users, particularly when simulating very complex individuals. I think that I will do this for the supporting information showing an advanced example.
Custom crop and pesticide initialisation and rotation
I have now introduced custom crop and pesticide initialisation and
rotation, as I set out to do in the 13 JUN
2022 notes below. The code appears to be stable, and I have
introduced two new automated tests in the package to make sure that the
feature remains stable. Below is an example of how to have a custom
rotation regime for crops and pesticides; I will start with crops.
First, I will use the mg_n1 gmatrix that comes with the
resevol package (I could have also made a dummy one; the
details are not important here).
data("gmatrices");Now I will run a simulation. The run_farm_sim takes some
simple options for crop_init and
crop_rotation_type. This includes
crop_init = "random", which randomly assigns each farm to
one of the potential crops. There are also three
crop_rotation_type arguments, including 1 (no
crop rotation), 2 (random crop choice whenever it is time
to change), and 3 (cycling through crops by number; i.e.,
1, 2, 3, 1, 2, 3). I might add more preset options like this later, but
here is how to run a simulation with two crops, random crop
initialisation on farms, and random rotation between crops.
sim <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4,
repro = "asexual", time_steps = 10, farms = 4,
crop_number = 2,
crop_init = "random",
crop_rotation_type = 2,
print_inds = FALSE, print_gens = FALSE,
print_last = FALSE, get_stats = FALSE);The new feature is that crop_init and
crop_rotation_type can take custom arguments, so
potentially any initialisation and cycling of crop choices is possible.
With crop_init, we can specify which farm gets what crop at
the start of the simulation using a fector of length farms.
If, for example, we want farms 1 and 4 to have crop 1, and farms 2 and 3
to have crop 2, we could run the below.
sim <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4,
repro = "asexual", time_steps = 10, farms = 4,
crop_number = 2,
crop_init = c(1, 2, 2, 1),
crop_rotation_type = 2,
print_inds = FALSE, print_gens = FALSE,
print_last = FALSE, get_stats = FALSE);Note that the position of farms on the landscape is shown below.
## [,1] [,2] [,3] [,4]
## [1,] 1 1 2 2
## [2,] 1 1 2 2
## [3,] 3 3 4 4
## [4,] 3 3 4 4Hence, this landscape would start out with farms numbered 1 and 4,
diagonal from one another, applying crop 1, and farms numbered 2 and 3,
also diagonal, applying crop 2. Now, suppose that we want users of crop
1 to switch to crop 2 with a probability of 0.75, and retain the use of
crop 1 with a probability of 0.25. Likewise, suppose that we want users
of crop 2 to switch to crop 1 with a probability of 0.75, and retain the
use of crop 2 with a probability of 0.25. We can now specify this with a
square matrix to input in place of the number for the argument
crop_rotation_type.
change_mat <- matrix(data = c(0.25, 0.75, 0.75, 0.25), nrow = 2, ncol = 2);## [,1] [,2]
## [1,] 0.25 0.75
## [2,] 0.75 0.25Note that matrix elements \(M_{i,j}\) specify the probability of
changing from crop \(i\) to crop \(j\), meaning that row elements must sum to
1. This can be binary too, as in the default
crop_rotation_type = 3.
## [,1] [,2]
## [1,] 0 1
## [2,] 1 0The matrix above cycles back and forth between crops 1 and 2.
Overall, this should give a lot more flexibility in terms of specifying
how crops cycle over time. We can input the change_mat
directly as an argument in run_farm_sim.
sim <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4,
repro = "asexual", time_steps = 10, farms = 4,
crop_number = 2,
crop_init = c(1, 2, 2, 1),
crop_rotation_type = change_mat,
print_inds = FALSE, print_gens = FALSE,
print_last = FALSE, get_stats = FALSE);Nevertheless, the cycling is restricted to Markov chains, meaning that it cannot include cycles of crops such as 1, 1, 2, 1, 1, 2, 1, 1, etc., where the cycling repeats complex patterns. This is mainly due to the awkwardness of trying to allow for such cycles given simulations that can have high numbers of time steps (i.e., I do not think it is worth figuring out a system that allows for these kinds of cycles as an option in the code).
Pesticide rotation options work exactly the same way as crop rotation options, so I will not spend to much time on them. Here is an example of running one.
sim <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4,
repro = "asexual", time_steps = 10, farms = 4,
pesticide_number = 2,
pesticide_init = c(1, 2, 2, 1),
pesticide_rotation_type = change_mat,
print_inds = FALSE, print_gens = FALSE,
print_last = FALSE, get_stats = FALSE);In conclusion, this new feature will allow for a very flexible range of different crop and pesticide cycling options.
Custom landscape option
I have also added an option to customise the farming landscape, also
as I set out to do in the 13 JUN 2022 notes
below. The customisation is on the top layer of the landscape that
differentiates independent farms. What is required is a matrix that
includes a sequence of natural numbers in all elements, so if there are
4 farms, then all matrix elements must be either 1, 2, 3, or 4. Beyond
this requirement, there is no restriction on where different farms are
placed; the do not even need to be contiguous on the landscape. Custom
terrains can be inserted using the terrain argument in
run_farm_sim. Below makes a custom terrain with three
farms.
custom_terrain <- matrix(data = 0, nrow = 8, ncol = 8);
custom_terrain[,1:2] <- 1;
custom_terrain[,2:6] <- 2;
custom_terrain[,7:8] <- 3;## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 2 2 2 2 2 3 3
## [2,] 1 2 2 2 2 2 3 3
## [3,] 1 2 2 2 2 2 3 3
## [4,] 1 2 2 2 2 2 3 3
## [5,] 1 2 2 2 2 2 3 3
## [6,] 1 2 2 2 2 2 3 3
## [7,] 1 2 2 2 2 2 3 3
## [8,] 1 2 2 2 2 2 3 3We can insert the terrain into the landscape with the
terrain argument.
sim <- run_farm_sim(mine_output = mg_n1, N = 100, xdim = 4, ydim = 4,
repro = "asexual", time_steps = 10, farms = 4,
terrain = custom_terrain,
print_inds = FALSE, print_gens = FALSE,
print_last = FALSE, get_stats = FALSE);Note that a custom terrain will override the arguments
farms, xdim, and ydim, so even
though the above sets farms = 4, there are actually only
three farms on the landscape. I am tempted to put up a warning for this,
but I am actually not sure if it is needed; my hope is that anyone
customising the landscape would recognise this from the
documentation.
Also note that these terrain values do not necessarily need to be farms. Through the use of a custom landscape and pesticide rotation option, these cells could represent something like diversionary feeding sites or even buildings or rivers. Consider, for example, the following matrix.
complex_terrain <- matrix(data = 0, nrow = 12, ncol = 12)
complex_terrain[1:2, 1:3] <- 1;
complex_terrain[3:7, 1:2] <- 1;
complex_terrain[8:12, 1:2] <- 2;
complex_terrain[1:2, 6:9] <- 3;
complex_terrain[3:6, 7:9] <- 3;
complex_terrain[1:6, 10:12] <- 4;
complex_terrain[10:12, 7] <- 5;
complex_terrain[7:12, 8:12] <- 5;
complex_terrain[3:12, 3:4] <- 6;
complex_terrain[7:9, 5] <- 6;
complex_terrain[1:2, 4:5] <- 7;
complex_terrain[3:6, 5:6] <- 7;
complex_terrain[7:9, 6:7] <- 7;
complex_terrain[10:12, 5:6] <- 7;## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1 1 1 7 7 3 3 3 3 4 4 4
## [2,] 1 1 1 7 7 3 3 3 3 4 4 4
## [3,] 1 1 6 6 7 7 3 3 3 4 4 4
## [4,] 1 1 6 6 7 7 3 3 3 4 4 4
## [5,] 1 1 6 6 7 7 3 3 3 4 4 4
## [6,] 1 1 6 6 7 7 3 3 3 4 4 4
## [7,] 1 1 6 6 6 7 7 5 5 5 5 5
## [8,] 2 2 6 6 6 7 7 5 5 5 5 5
## [9,] 2 2 6 6 6 7 7 5 5 5 5 5
## [10,] 2 2 6 6 7 7 5 5 5 5 5 5
## [11,] 2 2 6 6 7 7 5 5 5 5 5 5
## [12,] 2 2 6 6 7 7 5 5 5 5 5 5We can visualise the terrain a bit more clearly with an image below.
We can imagine the dark maroon colour as being a river (7), and the red region on the immediate left of the river as being a nature reserve (6). To get the simulations to run treating it this way, we would only need to spray a crop and pesticide that has no effect on the pests on the river (effectively making it an empty area of landscape in terms of feeding or resistance), and a pesticide but not crop in the nature reserve. First we can initialise the crop and pesticide vectors. Assume that there are 2 legitimate crops and pesticides that the actual farms (1-5) will use.
initial_crops <- c(1, 2, 1, 2, 1, 3, 4);
initial_pesticides <- c(1, 1, 2, 2, 1, 3, 3);We can then clarify the matrix so that there is only a cycle between 1 and 2, but never 3 and 4 for either the pesticide or the crops. First I will do the crops.
crop_cycle <- matrix(data = 0, nrow = 4, ncol = 4);
crop_cycle[1, 2] <- 1;
crop_cycle[2, 1] <- 1;
crop_cycle[3, 3] <- 1;
crop_cycle[4, 4] <- 1;## [,1] [,2] [,3] [,4]
## [1,] 0 1 0 0
## [2,] 1 0 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1Now I will do the pesticide.
pesticide_cycle <- matrix(data = 0, nrow = 3, ncol = 3);
pesticide_cycle[1, 2] <- 1;
pesticide_cycle[2, 1] <- 1;
pesticide_cycle[3, 3] <- 1;## [,1] [,2] [,3]
## [1,] 0 1 0
## [2,] 1 0 0
## [3,] 0 0 1So farms (1-5) will swap between crops 1 and 2, while the reserve (6) will maintain crop 3 and the river (7) will be ‘crop’ 4. And farms will swap between pesticides 1 and 2, while the reserve and river will hold pesticide 3. Then, to make crop 4 and pesticide 3 not really exist, we just make it impossible for the pests to consume these.
sim <- run_farm_sim(mine_output = mg_n1, N = 100,
repro = "asexual", time_steps = 10,
crop_number = 4, pesticide_number = 3,
terrain = complex_terrain,
food_consume = c(0.25, 0.25, 0.25, 0),
pesticide_consume = c(0.1, 0.1, 0),
crop_init = initial_crops,
crop_rotation_type = crop_cycle,
pesticide_init = initial_pesticides,
pesticide_rotation_type = pesticide_cycle,
print_inds = FALSE, print_gens = FALSE,
print_last = FALSE, get_stats = FALSE);The output shows that the pests leave the ‘food’ on the river alone (see layer 5), while the reserve (see layer 4) gets used very much.
## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1 1 1 7 7 3 3 3 3 4 4 4
## [2,] 1 1 1 7 7 3 3 3 3 4 4 4
## [3,] 1 1 6 6 7 7 3 3 3 4 4 4
## [4,] 1 1 6 6 7 7 3 3 3 4 4 4
## [5,] 1 1 6 6 7 7 3 3 3 4 4 4
## [6,] 1 1 6 6 7 7 3 3 3 4 4 4
## [7,] 1 1 6 6 6 7 7 5 5 5 5 5
## [8,] 2 2 6 6 6 7 7 5 5 5 5 5
## [9,] 2 2 6 6 6 7 7 5 5 5 5 5
## [10,] 2 2 6 6 7 7 5 5 5 5 5 5
## [11,] 2 2 6 6 7 7 5 5 5 5 5 5
## [12,] 2 2 6 6 7 7 5 5 5 5 5 5
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.75 0.75 0.25 0 0 0 1.00 1.00 1.00 0 0 0
## [2,] 1.00 0.00 0.00 0 0 0 0.75 1.00 1.00 0 0 0
## [3,] 0.75 0.25 0.00 0 0 0 0.75 1.00 1.00 0 0 0
## [4,] 0.75 0.00 0.00 0 0 0 1.00 1.00 1.00 0 0 0
## [5,] 0.75 0.00 0.00 0 0 0 0.25 1.00 0.75 0 0 0
## [6,] 0.75 0.00 0.00 0 0 0 0.00 1.00 1.00 0 0 0
## [7,] 0.50 0.25 0.00 0 0 0 0.00 1.00 1.00 1 1 1
## [8,] 0.00 0.00 0.00 0 0 0 0.00 0.25 0.25 1 1 1
## [9,] 0.00 0.00 0.00 0 0 0 0.00 0.50 1.00 1 1 1
## [10,] 0.00 0.00 0.00 0 0 0 0.00 0.25 1.00 1 1 1
## [11,] 0.00 0.00 0.00 0 0 0 0.75 1.00 1.00 1 1 1
## [12,] 0.00 0.00 0.00 0 0 0 0.00 0.75 1.00 1 1 1
##
## , , 3
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.00 0.00 0 0 0 0 0 0 0 1 1 1
## [2,] 0.00 0.00 0 0 0 0 0 0 0 1 1 1
## [3,] 0.00 0.00 0 0 0 0 0 0 0 1 1 1
## [4,] 0.00 0.00 0 0 0 0 0 0 0 1 1 1
## [5,] 0.00 0.00 0 0 0 0 0 0 0 1 1 1
## [6,] 0.00 0.00 0 0 0 0 0 0 0 1 1 1
## [7,] 0.00 0.00 0 0 0 0 0 0 0 0 0 0
## [8,] 1.00 0.00 0 0 0 0 0 0 0 0 0 0
## [9,] 0.25 0.00 0 0 0 0 0 0 0 0 0 0
## [10,] 0.25 0.00 0 0 0 0 0 0 0 0 0 0
## [11,] 0.75 0.25 0 0 0 0 0 0 0 0 0 0
## [12,] 0.50 1.00 0 0 0 0 0 0 0 0 0 0
##
## , , 4
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [9,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [10,] 0 0 0.0 0 0 0 0 0 0 0 0 0
## [11,] 0 0 0.5 0 0 0 0 0 0 0 0 0
## [12,] 0 0 0.0 0 0 0 0 0 0 0 0 0
##
## , , 5
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0 0 0 1 1 0 0 0 0 0 0 0
## [2,] 0 0 0 1 1 0 0 0 0 0 0 0
## [3,] 0 0 0 0 1 1 0 0 0 0 0 0
## [4,] 0 0 0 0 1 1 0 0 0 0 0 0
## [5,] 0 0 0 0 1 1 0 0 0 0 0 0
## [6,] 0 0 0 0 1 1 0 0 0 0 0 0
## [7,] 0 0 0 0 0 1 1 0 0 0 0 0
## [8,] 0 0 0 0 0 1 1 0 0 0 0 0
## [9,] 0 0 0 0 0 1 1 0 0 0 0 0
## [10,] 0 0 0 0 1 1 0 0 0 0 0 0
## [11,] 0 0 0 0 1 1 0 0 0 0 0 0
## [12,] 0 0 0 0 1 1 0 0 0 0 0 0Similarly, the pests will only be affected by pesticides 1 and 2, while 3 is effectively nothing. Note that this is of course a trick for modelling landscape areas without food or pesticide, but I do not see a reason to make this more explicit in the function (unless I receive requests due to confusion, but my hope is it is clear from these notes and documentation).
As I prepare the pre-print and software note for the resevol package, I am thinking that might be worth doing just a bit more coding so that the publication can include four more features:
- More explicit crop and pesticide rotation, as
mentioned in Issue #43. I
think that this would be fairly easy to implement and worth the time.
Currently, the
run_farm_simargumentcrop_rotation_typeaccepts either 1 or a 2. A 1 results in no rotation of crops on the landscape, and a 2 results in a random transition between types. I think that a custom transition matrix should be possible to insert, whereby each element should specify the probability of transitioning from crop i (row) to crop j (column). In the simplest case (e.g., could becrop_rotation_type = 3), this could just be a binary matrix that forms a cyclical graph. The current options 1 and 2 could even be subsumed in the transition matrix. - Allow custom landscapes. I think that this should
be possible by allowing users to put their own matrix in a function,
with numbers on the matrix representing different farms (integers
greater than zeros) or public land (zeros). A new
custom_landscapefunction could check this matrix, then use it as the top layer of an array. - Custom plotting functions. This one is tricky given the flexibility of simulations and therefore the complexity of the output, but it might be worth including at least some basic plotting functions, or summary functions.
- Simulated annealing algorithm for the
mine_gmatrixoption. This could be fairly easy to implement, and for some cases might actually find genome solutions faster than the evolutionary algorithm.
My hope is that these different tasks could be done fairly quickly, and that they will make writing up the paper go more smoothly.
There had for some time been an annoying issue of resevol CRAN checks having ERRORs on Windows systems with ix86 architecture. As it turns out, this was caused by some errors in specifying the memory for arrays, which was fixed in this commit. Just to make sure that everything is clear and correct, I have made the cast much clearer in a new commit. Tests using the rhub package and on Winbuilder confirm that everything works as intended. Valgrind found no memory leaks; everything appears fine, so I have updated the version to v0.2.0.9 and submitted the fix to CRAN.
I have now run the same analysis that I did on 19 AUG 2021, but with traits that evolve
independently. This was done by setting a covariance matrix to be the
identity matrix for four traits. Here is the mine_gmatrix
output that was used for the simulations.
## [[1]]
## [1] 12.00 6.00 2000.00 18000.00 0.01 0.01 3200.00 0.01
## [9] 1200.00 6.00 -9.00 0.10 0.00
##
## [[2]]
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1
##
## [[3]]
## [,1] [,2] [,3] [,4]
## [1,] 1.01323287 -0.91056130 0.22300841 0.978403744
## [2,] 0.27599802 -0.37879572 -1.01693690 0.762334090
## [3,] -0.03203048 -0.51511158 0.79037437 0.002639731
## [4,] -0.09193560 0.07907323 0.14790265 -0.039062606
## [5,] -0.73131494 0.07549138 0.17234648 -0.009426905
## [6,] -0.50890530 0.18561143 0.72008751 0.496708973
## [7,] -0.01427086 -0.28208810 -0.11196553 -0.537226514
## [8,] 0.25274376 0.15911871 -0.11647583 0.011376558
## [9,] -0.21570241 -0.33741651 0.47711231 -0.459624091
## [10,] -0.23138642 -0.92727329 0.04955606 -0.017246277
## [11,] -0.19430907 -0.38257361 -0.26531432 0.630049848
## [12,] -0.66903676 0.17674572 -0.16014889 0.067177474
##
## [[4]]
## , , 1
##
## [,1] [,2] [,3] [,4]
## [1,] 0.3666743 0.86477950 0.61013519 0.4585473
## [2,] -0.7287975 0.89967854 -0.26984724 0.1715412
## [3,] -0.2753835 -0.27174429 -0.02783652 0.5228947
## [4,] 0.2106197 -0.09700542 -0.95923363 0.2126882
##
## , , 2
##
## [,1] [,2] [,3] [,4]
## [1,] -0.5464374 -0.57925350 0.2394385 0.13922962
## [2,] -0.2060458 0.01186567 -0.2479210 0.52160749
## [3,] 0.6099824 -0.40519047 -0.1605294 -0.05806732
## [4,] 0.7619924 -0.03927812 1.0184070 0.76831379
##
## , , 3
##
## [,1] [,2] [,3] [,4]
## [1,] 0.01579773 -0.0176360 0.1042295 -0.6998017767
## [2,] 1.12913619 -0.1744167 0.1694019 -0.0537070606
## [3,] 0.11367376 -0.5660887 0.5790242 -0.0007441308
## [4,] 0.00766533 -0.7063434 -0.6071945 -0.0742630499
##
## , , 4
##
## [,1] [,2] [,3] [,4]
## [1,] -0.4415180 -0.45882179 0.189424993 -0.4583560
## [2,] 0.4881191 -0.07763796 0.534154371 -0.3413463
## [3,] -0.1194739 0.71146855 0.622996667 0.7601848
## [4,] -0.2681222 0.75851682 0.004319161 -0.7493150
##
## , , 5
##
## [,1] [,2] [,3] [,4]
## [1,] 0.3594965 -0.44704612 0.1057529 0.6412851
## [2,] -0.4734836 -0.02154132 0.8362356 0.0703729
## [3,] -0.1015575 0.59180270 -0.2234910 0.3626638
## [4,] 0.7883569 0.48251087 0.6203973 0.1364310
##
## , , 6
##
## [,1] [,2] [,3] [,4]
## [1,] 0.6413329 -0.0520821 0.0313601 -0.7862978
## [2,] 0.2766572 0.5259953 1.0408831 0.7330976
## [3,] -0.4612240 0.6870955 -0.2413829 -0.1735197
## [4,] -0.8082452 -0.7595234 0.7075496 -0.6229779
##
##
## [[5]]
## [,1] [,2] [,3] [,4]
## [1,] -0.11900916 0.72971816 -0.33534413 -0.16856932
## [2,] -0.56655781 -0.01421153 -0.44314434 0.25331983
## [3,] 0.02974235 0.39626672 0.35083549 -0.16975625
## [4,] 0.08175985 0.01612911 0.06762713 0.03284772
## [5,] -0.10161875 -0.04482800 0.34340020 0.33128165
## [6,] 0.27189293 0.34943436 0.27667359 0.51511528
## [7,] -0.22418646 -0.17966680 0.14259016 -0.36500830
## [8,] 0.11894101 -0.03845361 -0.18004009 -0.05569821
## [9,] -0.06241772 0.05742287 0.38727606 -0.27966268
## [10,] -0.59027383 0.15893478 0.33824411 -0.16558486
## [11,] -0.38817067 0.16736336 -0.01299629 0.34651855
## [12,] -0.16592813 -0.15620554 0.18829157 0.39244612
##
## [[6]]
## [,1] [,2] [,3] [,4]
## [1,] 1.06913329 -0.072243449 0.064732772 -0.056487174
## [2,] -0.07224345 0.874066557 0.007063041 -0.007362562
## [3,] 0.06473277 0.007063041 0.948800593 -0.009399430
## [4,] -0.05648717 -0.007362562 -0.009399430 1.001389459
##
## [[7]]
## [1] 1.0132328710 -0.9105612976 0.2230084069 0.9784037440 0.2759980196
## [6] -0.3787957167 -1.0169369015 0.7623340902 -0.0320304800 -0.5151115759
## [11] 0.7903743678 0.0026397307 -0.0919356017 0.0790732281 0.1479026489
## [16] -0.0390626056 -0.7313149438 0.0754913824 0.1723464850 -0.0094269050
## [21] -0.5089053016 0.1856114264 0.7200875059 0.4967089726 -0.0142708566
## [26] -0.2820881047 -0.1119655280 -0.5372265136 0.2527437636 0.1591187075
## [31] -0.1164758339 0.0113765575 -0.2157024076 -0.3374165063 0.4771123077
## [36] -0.4596240914 -0.2313864218 -0.9272732890 0.0495560606 -0.0172462772
## [41] -0.1943090688 -0.3825736088 -0.2653143217 0.6300498485 -0.6690367591
## [46] 0.1767457169 -0.1601488935 0.0671774739 0.3666742930 0.8647795015
## [51] 0.6101351924 0.4585473120 -0.7287975411 0.8996785444 -0.2698472400
## [56] 0.1715411534 -0.2753834548 -0.2717442880 -0.0278365186 0.5228947384
## [61] 0.2106197266 -0.0970054220 -0.9592336331 0.2126882478 -0.5464374357
## [66] -0.5792534966 0.2394384985 0.1392296218 -0.2060458352 0.0118656721
## [71] -0.2479210148 0.5216074862 0.6099824149 -0.4051904672 -0.1605293560
## [76] -0.0580673216 0.7619923905 -0.0392781166 1.0184070112 0.7683137886
## [81] 0.0157977321 -0.0176359999 0.1042294944 -0.6998017767 1.1291361885
## [86] -0.1744166518 0.1694018671 -0.0537070606 0.1136737596 -0.5660887218
## [91] 0.5790242132 -0.0007441308 0.0076653298 -0.7063434037 -0.6071945083
## [96] -0.0742630499 -0.4415179511 -0.4588217919 0.1894249926 -0.4583559714
## [101] 0.4881190547 -0.0776379596 0.5341543712 -0.3413463261 -0.1194739358
## [106] 0.7114685497 0.6229966674 0.7601847997 -0.2681222285 0.7585168195
## [111] 0.0043191611 -0.7493149802 0.3594964704 -0.4470461217 0.1057529089
## [116] 0.6412850683 -0.4734836453 -0.0215413243 0.8362356329 0.0703728975
## [121] -0.1015575348 0.5918026977 -0.2234910091 0.3626638219 0.7883569250
## [126] 0.4825108655 0.6203972867 0.1364309830 0.6413329012 -0.0520821005
## [131] 0.0313601037 -0.7862978038 0.2766571825 0.5259952870 1.0408830938
## [136] 0.7330976036 -0.4612239907 0.6870955302 -0.2413828743 -0.1735197274
## [141] -0.8082452080 -0.7595233690 0.7075495824 -0.6229778721
##
## [[8]]
## [1] -5.625305Note that the parameter values for mining the matrix were identical to those chosen in the simulations from 19 AUG 2021. The results are as expected; all traits evolve in response to selection to feed on both crop types and resist both pesticide types.
Asexual, no rotation
