bug in simulation function?

cecileane commented 4 years ago

I am afraid there is a bug with the simulation function, when nreps is more than 1. Example below, with a 3-taxon tree, 2 independent traits, both with very small within-species variation to see the problem.

> tree = read.tree(text="((A:0.4,B:0.4):4.4,C:5.2);")
> trait_cov = matrix(c(1,0, 0,2), 2,2)
> trait_cov
     [,1] [,2]
[1,]    1    0
[2,]    0    2
> within_species_cov = diag(c(1e-10,1e-10))
> within_species_cov
      [,1]  [,2]
[1,] 1e-10 0e+00
[2,] 0e+00 1e-10

Now to the simulation:

> set.seed(1) # for reproducibility
> sim_data_intra <- simtraits(v=trait_cov,
+                             intraspecific=within_species_cov,
+                             tree=tree, nreps=2, nmissing=0)
> sim_data_intra$original_X
          V1        V2
A -0.2308388 -2.039192
B  0.4964700 -1.314619
C -0.6836394  6.144626

It looks good so far: V1 and V2 have "reasonable" simulated species means. The problem comes next. The individual values are good for V1: each individual has a value almost equal to the species mean (up the the 5th decimal point, because of within-species variance of 1e-10). But the second variable V2 seems to be messed up:

for the 1st individual of each species, V2 ≈ V1 --that's the issue.
for the 2nd individual of each species, V2 seems okay: equal to the species mean of V2 (from above) up to the 5th decimal point, roughly.

> sim_data_intra$trait_data
  species         V1         V2
1       A -0.2308539 -0.2308427
2       B  0.4964587  0.4964704
3       C -0.6836476 -0.6836453
4       A -0.2308427 -2.0391954
5       B  0.4964704 -1.3146416
6       C -0.6836453  6.1446354

ericgoolsby commented 4 years ago

I believe you're correct! Line 92 of phylopars_main.R (version 0.2.11) reads:

X[[j]][1:ntaxa + (jj-1)(ntaxa),] <- simdat_j[,1:ntraitsjj]

but it should be

X[[j]][1:ntaxa + (jj-1)(ntaxa),] <- simdat_j[,(1:ntraits-1)nreps+jj]

I believe this will correct the problem. For the example above, the revised function gives:

  species         V1        V2
1       A -0.2308539 -2.039187
2       B  0.4964587 -1.314626
3       C -0.6836476  6.144626
4       A -0.2308427 -2.039195
5       B  0.4964704 -1.314642
6       C -0.6836453  6.144635

This will be included in version >=0.2.12. Thanks for catching this!

cecileane commented 4 years ago

Thank you for your diagnosis and fix, Eric!

ericgoolsby / Rphylopars

bug in simulation function? #33