Notes

A PDF version is provided, but this is not a very useful format.  In the Mathematica notebook, ancillary runs are hidden (double-click brackets on the right to see them).

First, Evaluate Initialisation Cells. The code needed to generate the figures is provided nearby,  but some simulation output may need to be read in from auxiliary files to reproduce the exact output:

“Pedigrees N=30 neutral II”
“Initial population N=30 II”
“Pairwise F N=30 neutral II”
“p, zA, zD, N=30 neutral II”
“p, zA, zD, N=30 neutral III”
“p, zA, zD, N=30 neutral IV”
"p, zA, zD, N=30 neutral III";
"p, zA, zD, N=30 selected k=5";
“Neutral identities 16th Jan”
“Pairwise F N=30 neutral II”
“selected reps N=30 selected k=1,2 300 reps III”

Some code (immediately above figures) is hidden: Use “Cell Properties>Open” to reveal.

The first  section summarises how to calculate identity coefficients from a pedigree; see ‘Algorithm …” for a fuller explanation.

Identity coefficients

Summary notes

Even with a single locus, the general algorithm for identities amongst genes in diploid individuals is non-trivial.  There are papers on this, but they give examples for only small pedigrees, and may assume that founders are equivalent.  Here, I describe the algorithm and attempt to implement it.

Genes are defined by individual and gene: 𝕚={i,1} or {i,2}, say.  We work with the probabilities that certain sets of genes are IBD. For example, F[{{(Subscript[𝕚, 1])},{(Subscript[𝕚, 2]),(Subscript[𝕚, 3])},{(Subscript[𝕚, 4]),(Subscript[𝕚, 5])}}] is the probability that these 3 sets of genes each trace back to 3 distinct founders.

Note that singletons are informative: F[{{(Subscript[1, 1]),(Subscript[1, 2])}}]=F[{{(Subscript[1, 1]),(Subscript[1, 2])},{(Subscript[2, 1])}}]+F[{{(Subscript[1, 1]),(Subscript[1, 2]),(Subscript[2, 1])}}]

Given some pedigree, the identity (Subscript[F, t]) is a linear combination of the (Subscript[F, t-1]), in which configurations change as a result of shared parents.
	- a single gene traces back to one or other parent with equal probability; two genes in the same individual must trace back to the two  parents, though those may be the same individual if there is selfing.  This generates a sum over possible sets of parents.
	- Configurations may now change as a result of shared parentage.  It should be that the order in which we consider shared parents makes no difference.
	- For each shared parent in turn, assign descendant lineages to one or other of two sets, with equal probability, and merge sets accordingly.
	- This gives a recursion back to the founder generation, when all identity coefficients are zero.

It should be possible to calculate all necessary probabilities by combining these identity coefficients.

The algorithm:
	- Make sure the list is ordered, at two levels.
	- at time zero, set to 1 if the set is empty or contains only singletons, zero otherwise
	- tracing backwards, set up all possible assignments of parents to the sets of genes, using genesToParents
	- for each parent individual, identify the sets from which incoming lineages come
	- assign these to all possible subsets of two (generating a potentially long list of possibilities)
	- if either of these subsets includes a mix of sets, set to zero
	- if incoming lineages all come from the same set, assign them to gene 1 or gene 2.

References

Abney M., 2009 A graphical algorithm for fast computation of identity coefficients and generalized kinship coefficients. Bioinformatics 25: 1561–1563.
	Not entirely clear, but described for two individuals only. 
Abney M., McPeek M. S., Ober C., 2000 Estimation of Variance Components of Quantitative Traits in Inbred Populations. The American Journal of Human Genetics 66: 629–650.
	discusses how one may estimate the various dominance variance components. Based on Karigl’s algorithm for the 9 condensed identity coefficients, but implemented more efficiently; includes a short appendix C on the CLT with dominance.
García-Cortés L. A., 2015 A novel recursive algorithm for the calculation of the detailed identity coefficients. Genetics selection evolution 47: 1319–4.
Karigl G., 1982 A mathematical approach to multiple genetic relationships. Theor. Popul. Biol.: 1–15.
Karigl G., 1981 A recursive algorithm for the calculation of identity coefficients. Annals Of Human Genetics 45: 290–305.
Kirkpatrick B., Ge S., Wang L., 2018 Efficient computation of the kinship coefficients (O Stegle, Ed.). Bioinformatics 35: 1002–1008.
	Efficient algorithm, suitable for large pedigrees; described only for two individuals.
Thompson E. A., 2013 Identity by Descent: Variation in Meiosis, Across Genomes, and in Populations. Genetics 194: 301–326.
	Overview of the relation between IBD and the ARG

An example

This sets up a random WF pedigree (no selfing)

pedigreeList[𝒫]=Table[Parents[randomPedigree[{10,10}]],{50}];

The mean IBD between the two genes within an individual increases at about the right rate. (Actually, because I do  not allow selfing, this should increase IBD above 20).

tt1=Table[Mean[Fx[𝒫,t][{{{#,1},{#,2}}}]&/@Range[10]],{t,0,50}]//N

-1/D[Fit[Drop[Log[1-tt1],20],{1,t},t],t]

Show[ListLogPlot[1-tt1],LogPlot[(((1-(1/20)))^t),{t,0,50},PlotStyle->Red]]

tt2=Table[Mean[Flatten[Table[Fx[𝒫,t][{{{i,1},{j,1}}}],{i,2,10},{j,1,i-1}],1]],{t,0,50}]

-1/D[Fit[Drop[Log[1-tt2],20],{1,t},t],t]

Show[ListLogPlot[1-tt2],LogPlot[(((1-(1/20)))^t),{t,0,50},PlotStyle->Red]]

This is the chance of three-way IBD between two genes in one diploid and another gene in another diploid, at times 10, 20:

Fx[𝒫,10][{{{1,1},{1,2},{2,1}}}]

Fx[𝒫,20][{{{1,1},{1,2},{2,1}}}]

Simulations for the paper

Setting up the pedigree and initial population

A single pedigree, and a single initial population of genotypes, is chosen, labelled 𝒩.  Then, replicates will be evolved, possibly with within-family selection but on a fixed pedigree.  Neutral simulations will be denoted 𝒩[j] for replicate j, and simulations that select the best of k within families will be denoted 𝒩[j]

Note that storing genotypes uses a lot of memory and may not be feasible for all replicates.  Storing multi-way identity coefficients may also run into difficulties.

Setting up the initial population: N=30, M=1000 loci, γ=2

This sets up a WF population (but with no selfing) of N=30 diploid individuals, and M=1000 loci. Allelic effects are all α=(1/(Sqrt[M])), and dominance is h=0,1; probability of h=1 is (1/2)+(2/(Sqrt[M])), which should ensure (Subscript[V, D]) of order 1.  Allele frequencies are sampled from a large base population with allele frequencies distributed as (1/(p q)).  The reference frequencies are those in the actual sampled population.

ninds=30;nloci=1000;
αh=setαh[nloci,2.]//N;
jr=randomJ[2ninds,0.5,0.999,nloci];pr=jr/(2ninds);zb=zbar[pr,αh];
pop[𝒩]={initPop[ninds,jr]};

Note that with this model, the maximum range of the additive component is ±(Sqrt[(Subscript[n, loci])])=±31.

Save["Pedigrees N=30 neutral",{pedigree,pedigreeList}];

initialP=pop[𝒩];(* avoid storing all the simulated genotypes *)
Save["Initial population N=30",{initialP,αh}];

Making the pedigrees

This makes the pedigrees, over T=100 generations. There is no selfing:

T=100;
makePedigrees[𝒩,ninds,T];

This calculates pairwise coefficients within and between individuals:

pfw[t_]=Table[Fx[𝒩,t][{{{i,1},{i,2}}}],{i,ninds}];
pfb[t_]:=Flatten[Table[Fx[𝒩,t][{{{i,1},{j,1}}}],{i,2,ninds},{j,1,i-1}]];

The identity within individuals lags that between by one generation (black, blue).  The mean identity between genes in different individuals follows 1-(((1-(1/(2N))))^t) well (red). (I ignore the slight adjustment due to no selfing):

tt=Table[{t,Mean[pfb[t]]},{t,1,50}];
ttw=Table[{t,Mean[pfw[t]]},{t,1,50}];
Show[ListLinePlot[tt],ListLinePlot[ttw,PlotStyle->Black],Plot[1-(((1-(1/(2ninds))))^t),{t,0,50},PlotStyle->Red]]

This shows histograms of identity between individuals at t=10, 25, 50 generations.

cf= ColorData["DarkRainbow"];
Show[Histogram[pfb[#1],20,ChartStyle->cf[#/50],PlotRange->{{0,0.7},{0,100}}]&/@{10,25,50}]

Neutral, with 100 replicates

A single pedigree, and a single initial population of genotypes, is chosen, labelled 𝒮.  Then, replicates will be evolved, possibly with within-family selection but on a fixed pedigree.  Neutral simulations will be denoted 𝒩[j] for replicate j, and simulations that select the best of k within families will be denoted 𝒮[j]

Note that storing genotypes uses a lot of memory and may not be feasible for all replicates.  Storing multi-way identity coefficients may also run into difficulties.

Setting up the initial population: N=30, M=1000 loci, γ=2

This sets up a WF population (but with no selfing) of N=30 diploid individuals, and M=1000 loci. Allelic effects are all α=(1/(Sqrt[M])), and dominance is h=0,1; probability of h=1 is (1/2)+(2/(Sqrt[M])), which should ensure (Subscript[V, D]) of order 1.  Allele frequencies are sampled from a large base population with allele frequencies distributed as (1/(p q)).  The reference frequencies are those in the actual sampled population.

ninds=30;nloci=1000;
αh=setαh[nloci,2.]//N;
jr=randomJ[2ninds,0.5,0.999,nloci];pr=jr/(2ninds);zb=zbar[pr,αh];
pop[𝒩]={initPop[ninds,jr]};

Note that with this model, the maximum range of the additive component is ±(Sqrt[(Subscript[n, loci])])=±31.

Making the pedigrees

This makes the pedigrees, over T=100 generations. There is no selfing:

T=50;
makePedigrees[𝒩,ninds,T];

This calculates pairwise coefficients within and between individuals:

pfw[t_]=Table[Fx[𝒩,t][{{{i,1},{i,2}}}],{i,ninds}];
pfb[t_]:=Flatten[Table[Fx[𝒩,t][{{{i,1},{j,1}}}],{i,2,ninds},{j,1,i-1}]];

The identity within individuals lags that between by one generation (black, blue).  The mean identity between genes in different individuals follows 1-(((1-(1/(2N))))^t) well (red). (I ignore the slight adjustment due to no selfing):

tt=Table[{t,Mean[pfb[t]]},{t,1,50}];
ttw=Table[{t,Mean[pfw[t]]},{t,1,50}];
Show[ListLinePlot[tt],ListLinePlot[ttw,PlotStyle->Black],Plot[1-(((1-(1/(2ninds))))^t),{t,0,50},PlotStyle->Red]]

This shows histograms of identity between individuals at t=10, 25, 50 generations.

cf= ColorData["DarkRainbow"];
Show[Histogram[pfb[#1],20,ChartStyle->cf[#/50],PlotRange->{{0,0.7},{0,100}}]&/@{10,25,50}]

Saving pedigrees and initial population

Save["Pedigrees N=30 neutral II",{pedigree,pedigreeList}];

initialP=pop[𝒩];(* avoid storing all the simulated genotypes *)
Save["Initial population N=30 II",{initialP,αh}];

Save["Pairwise F N=30 neutral II",Fx];

Initial values of (Subscript[V, a]),(Subscript[V, d]),cov[A,D],ℐ, (Subscript[σ, DI]^2), (Subscript[σ, ADI])

These are the initial components, measured in the reference population.Note that (SuperStar[ι])=(Subscript[V, D]) for biallelic loci.

{Va[pr,αh],Vd[pr,αh],Vad[pr,αh],ℐ[pr,αh],ℐ2[pr,αh],σADI[pr,αh],σ2DI[pr,αh] }

There is a subtle issue here, however: these will not be precisely the same as for genes sampled from the single actual base population. For example, 𝔼[(SuperStar[ζζ])] equals p q for identical genes, but equals -((p q)/(2N-1)) if two distinct genes are sampled from a finite population with known frequency p.

Evolving 1000 neutral replicates, conditional on the pedigree

Setting up

This sets up the replicates - and is very fast.  However, the population genotypes take up a lot of memory, which limits the # of replicates that can be safely generated.

Dynamic[r]

nreps = 1000;
Timing[Do[evolve[{𝒩, 𝒩[r]}];
   pList[𝒩[r]] = N[Mean[Flatten[#, 1]]] & /@ pop[𝒩[r]];
   zAlist[𝒩[r]] = zA[#, pr, αh] & /@ pop[𝒩[r]];
   zDlist[𝒩[r]] = zD[#, pr, αh] & /@ pop[𝒩[r]];
   pop[𝒩[r]] =., {r, nreps}];]

{8699.35, Null}

Save["p, zA, zD, N=30 neutral II", {pList, zAlist, zDlist}];

Distributions of A, D in the final generation

These are the distributions of A, D in the final generation:

al=Table[Mean[zAlist[𝒩[r]]//Last],{r,1000}];
dl=Table[Mean[zDlist[𝒩[r]]//Last],{r,1000}];
GraphicsRow[Histogram/@{al,dl}]
{Mean[#],Variance[#]}&/@{al,dl}//TableForm

Mean (OverBar[A]), (OverBar[D]) through time

This shows the mean (OverBar[A]),(OverBar[D]) through time, averaged over 1000 replicates. (OverBar[D]) is compared to ι (Subscript[(OverBar[F]), w]) (red).  The mean (OverBar[A]) fluctuates close to zero, as expected.  The mean (OverBar[D]) decreases as the population becomes more inbred.  There are two discrepancies. First, at t=0 only, there is a strong negative (OverBar[D]), due to a chance deviation from HW in the single initial population.  Second, there is a constant discrepancy which arises because 𝔼[ζ (SuperStar[ζ])]=-p q/(2N-1) between two distinct genes that are sampled from the initial population. The lower graph shows that this discrepancy is as expected except in the last few generations.

GraphicsRow[{ListLinePlot[Transpose[{Range[0,50],Mean[Table[Mean/@zAlist[𝒩[r]],{r,1000}]]}],AxesLabel->{"t","(OverBar[A])"}],Show[ListLinePlot[Transpose[{Range[0,50],Mean[Table[Mean/@zDlist[𝒩[r]],{r,1000}]]}]],
ListLinePlot[Table[{t,ℐ[pr,αh] Mean[pfw[t]]},{t,0,50}],PlotStyle->Red],AxesLabel->{"t","(OverBar[D])"},PlotRange->All]}]

Show[ListLinePlot[Mean[Table[Mean/@zDlist[𝒩[r]],{r,1000}]]-Table[ℐ[pr,αh] Mean[pfw[t]],{t,0,50}]],Plot[-ℐ[pr,αh]/59.,{t,0,50},PlotStyle->Red]]

A, D vs F

The relation between (Subscript[F, 12]) and A,D (left, right), in generations 5, 10, 20, 40 (black, blue, purple, red).  SInce the same pedigree is used for all replicates, each individual is characterised by a single (Subscript[F, 12]).  The plots show A, D, averaged over 1000 replicates, for each individual; thus, there are N=30 points in each generation.  There is no systematic effect of F on the mean A, whereas the mean D decreases as ℐ F = -0.53F (red line at right).

gr[f_,l_]:=Show[MapThread[ListPlot[Transpose[{pfw[#1],Mean[Table[f[𝒩[r]][[#1+1]],{r,nreps}]]}],PlotStyle->#2]&,{{5,10,20,40},{Black,Blue,Purple,Red}}],PlotRange->All,AxesLabel->{"(Subscript[F, 12])",l}];
GraphicsRow[{gr[zAlist,"A"],Show[gr[zDlist,"D"],Plot[ℐ[pr,αh] f,{f,0,0.55},PlotStyle->Red]]}]

Change in mean and variance

All components are measured relative to the base population. Top row: means across 1000 replicates; bottom row: 10 replicates, individually. Left: mean (black), additive component (grey), dominance component (blue); right: actual variance (including LD); (Subscript[V, A]) (grey); (Subscript[V, D]) (blue); cov(A,D) (purple); (Subscript[V, A])+(Subscript[V, D])+2cov[A,D] (red).   Examining individual replicates (bottom row) we see that indeed, the actual (Subscript[V, G])(including LD) is much more variable than the underlying genic components.  Note that because I discard genotype frequencies, I cannot reconstruct the genic components of variance (which would be much smoother).

zbarList[r_]:=Mean/@(zAlist[𝒩[r]]+zDlist[𝒩[r]]);
varz[r_]:=Variance/@(zAlist[𝒩[r]]+zDlist[𝒩[r]]);
zal[r_]:=Mean/@zAlist[𝒩[r]];
zdl[r_]:=Mean/@zDlist[𝒩[r]];
val[r_]:=Variance/@zAlist[𝒩[r]];
vdl[r_]:=Variance/@zDlist[𝒩[r]];
vadl[r_]:=MapThread[Covariance,{zAlist[𝒩[r]],zDlist[𝒩[r]]}];
GraphicsGrid[{{
Show[MapThread[ListLinePlot[#1,PlotStyle->#2,PlotRange->{{0,T+1},{-1.5,1}},AxesLabel->{"time","(OverBar[z])-(Subscript[(OverBar[z]), 0]),(OverBar[A]),(OverBar[D])"}]&,
{Mean[Table[{zbarList[r],zal[r],zdl[r]},{r,nreps}]],{Black,Gray,Blue}}]],Show[MapThread[ListLinePlot[#1,PlotStyle->#2,PlotRange->{{0,T+1},{-0.1,0.5}},AxesLabel->{"time","V,(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])"}]&,{Mean[Table[{varz[r],val[r],vdl[r],val[r]+vdl[r]+2vadl[r],vadl[r]},{r,nreps}]],{Black,Gray,Blue,Red,Purple}}]]},
{Show[Table[MapThread[ListLinePlot[#1,PlotStyle->#2,PlotRange->{{0,T+1},{-1.5,1}},AxesLabel->{"time","(OverBar[z])-(Subscript[(OverBar[z]), 0]),(OverBar[A]),(OverBar[D])"}]&,
{{zbarList[r],zal[r],zdl[r]},{Black,Gray,Blue}}],{r,10}]],Show[Table[MapThread[ListLinePlot[#1,PlotStyle->#2,PlotRange->{{0,T+1},{-0.1,0.5}},AxesLabel->{"time","V,(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])"}]&,{{varz[r],val[r],vdl[r],val[r]+vdl[r]+2vadl[r],vadl[r]},{Black,Gray,Blue,Red,Purple}}],{r,10}]]}}]

Another attempt, starting with different replicated initial populations

Figures marked * to be in the paper

Evolving 1000 replicates, conditional on the pedigree

There was an obscure problem with running more than ~300 replicates - before, I managed 1000 with no problem.  The trick is, not to use Dynamic[r]...

<<"Pedigrees N=30 neutral II";
<<"Initial population N=30 II";

pr=Mean[Flatten[initialP[[1]],1]]//N;
nreps =1000;ninds=30;nloci=1000;T=50;
Timing[Do[
jrr=RandomInteger[BinomialDistribution[2ninds,#]]&/@pr;
pop[𝒩]={initPop[ninds,jrr]};
evolve[{𝒩,𝒩[r]}]; 
 pList[𝒩[r]] = N[Mean[Flatten[#, 1]]] & /@ pop[𝒩[r]];
zAlist[𝒩[r]] = zA[#, pr, αh] & /@ pop[𝒩[r]];
zDlist[𝒩[r]] = zD[#, pr, αh] & /@ pop[𝒩[r]];
   pop[𝒩[r]] =., {r, nreps}];]

Save["p, zA, zD, N=30 neutral III", {pList, zAlist, zDlist}];(* 300 replicates *)

Save["p, zA, zD, N=30 neutral IV", {pList, zAlist, zDlist}];(* 1000 replicates *)

Recovering the above…

(* sets up Pedigree, PedgreeList, initialP, αh, pList, zAlist, zDlist *)
<<"Pedigrees N=30 neutral II";
<<"Initial population N=30 II";
<<"p, zA, zD, N=30 neutral IV";

Distributions of A, D in the final generation

These are the distributions of A, D in the final generation:

al=Table[Mean[zAlist[𝒩[r]]//Last],{r,1000}];
dl=Table[Mean[zDlist[𝒩[r]]//Last],{r,1000}];
GraphicsRow[Histogram/@{al,dl}]
{Mean[#],Variance[#]}&/@{al,dl}//TableForm

Checking for an upward drift - not significant

tt=Table[zAlist[𝒩[r]][[40+1]],{r,nreps}];xx=Mean[tt];
Sort[xx]

There seems to be a significant upward drift of A, with 7/30<0 at t=40.

td=((Mean[xx])/(Sqrt[Variance[xx]/ninds]));
st=2 Sum[PDF[BinomialDistribution[30,1/2],j],{j,0,7}]//N;
{td,CDF[StudentTDistribution[29],-td],st}

However, the 30 points are correlated, and so the proper test is rather, whether the mean across all 30 is significant. The t-statistic for this is only 0.74, so here is no significant drift.  Note that the sd of (OverBar[A]) across replicates is ~0.5, much greater than the mean A for an one individual.

yy=Mean/@tt;
{Mean[yy],Variance[yy],((Mean[yy])/(Sqrt[Variance[yy]/1000]))}

ListPlot[Table[{i,#}&/@tt[[All,i]],{i,30}]]

** Predicting var(A), var(D), cov(A,D) across individuals

The IM predicts a multivariate normal distribution of the individuals in the pedigree, with mean and covariance that depend on the various identity coefficients.  Because we run replicates across a fixed pedigree, sampling from a given reference population. these can be predicted from the theory.  One can present this, as below, by plotting the variances against (say) (Subscript[F, i,j]).  I expect that there will be a systematic error, because as the population becomes more inbred, all the coefficients increase, and so a regression on one of the coefficients will be steeper than expected from just that coefficient. The noise around the regression will be due to random variation in the other F, plus the inevitable sampling error arising from estimation of a variance from 1000 values.

A more rigorous approach would be to plot the simulated variance against the prediction, based on all the F’s.  The variance around this should be a (χ^2) distribution with ~1000 d.f.. One can then ask how much of the variation is accounted for by pairwise F, and how much better the prediction is when all the F are included.  This should be feasible, though calculating fourth order identities for all generations may be challenging: there are 27405 4 way combinations to store for each of the 50 generations, so I guess here will be ~ 2 million coefficients altogether.

I simulate 50 generations each with 30 individuals, and so one can test 1500 variances and 1124250 covariances.  That would be excessive, but one might at least check covariances within chosen generations.

One could also check the shape of the distribution - but I propose checking normality only for the families generated below.

This constructs tables of the variance & covariance of A, D across replicates, for each of the 30 individuals, in a given generation:

** Predictions for families

In principle, the above provides a full check on the infinitesimal with dominance. However, it is easier to understand & present by focussing on families, generated from each pair in (say) generation 20.  The results can then be compared with the equations at the beginning of the paper.  I would also include a check on normality of the distributions of both shared and residual components.

The predictions are:

Needs sorting out!!

Residuals:
	var[A]=(1-(1/2)((Subscript[F, 11])+(Subscript[F, 22])))(Subscript[V, A])
cov[A,D]=((Subscript[F, 12])-(1/2)((Subscript[F, 112])+(Subscript[F, 122])))(Subscript[σ, ADI])
var[D]=(1/4)(3-(Subscript[F, 11])-(Subscript[F, 22])-(Subscript[F, 1122])-(Subscript[(OverTilde[F]), 1122])-(Subscript[(OverTilde[F]), 1212]))(Subscript[V, D])+(1/4)(3(Subscript[F, 12])-(Subscript[F, 1122])-(Subscript[F, 112])-(Subscript[F, 122]))(Subscript[σ, DI]^2)

Shared:
	𝔼[A]=0      𝔼[D]=ι (Subscript[F, 12])
	var[A]=(1/2)(1+(1/2)((Subscript[F, 11])+(Subscript[F, 22]))+2(Subscript[F, 12]))(Subscript[V, A])
	cov[A,D]=((Subscript[F, 12])+(1/2)((Subscript[F, 112])+(Subscript[F, 122])))(Subscript[σ, ADI])
	var[D] =(1/4)(1 + (Subscript[F, 11]) + (Subscript[F, 22]) + (Subscript[F, 1122])+  (Subscript[(OverTilde[F]), 1122]) +(Subscript[(OverTilde[F]), 1212])-4(Subscript[F, 12]^2)) (Subscript[V, D])+(1/4)((Subscript[F, 12])+(Subscript[F, 112])+(Subscript[F, 122])+ (Subscript[F, 1122]))(Subscript[σ, DI]^2)

Sum:
	var[A]=(1+(Subscript[F, 11]))(Subscript[V, A])
	cov[A,D]=2(Subscript[F, 12])(Subscript[σ, ADI])
	var[D] = (1 - (Subscript[F, 12]^2)) (Subscript[V, D])+ (Subscript[F, 12])(Subscript[σ, DI]^2)

Making and storing families of 1000 offspring from 100 random pairs in ten replicates, at generation 20

I will need several replicates, but do not need to make all 465 possible families within replicates.  Making 100 families would typically give ~6 selfed families.

pars=RandomSample[Flatten[Table[{j,i},{i,ninds},{j,i}],1],100]

This generates families of 1000 offspring, and stores {A,D} for them.

pedigree[𝒩20]=Take[pedigree[𝒩],20];
pop[𝒩20]=pop[𝒩];
makeF[r_][{i_,j_}]:=Module[{fm},
fm=makeFamily[pop[ℱ[r]][[21,i]],pop[ℱ[r]][[21,j]],1000];
Transpose[{zA[fm,pr,αh],zD[fm,pr,αh]}]];
Do[evolve[{𝒩20,ℱ[r]}];
zAlist[ℱ[r]] = zA[#, pr, αh] & /@ pop[ℱ[r]];
zDlist[ℱ[r]] = zD[#, pr, αh] & /@ pop[ℱ[r]];
fams[r]=makeF[r]/@pars;
pop[ℱ[r]]=.;,{r,10}];

Save["Neutral families 16th Jan",{fams,pars,zAlist,zDlist}];

It is quick to set up tables of family means and midparent values:

Do[cvf[r]=Covariance/@fams[r];
mnf[r]=Mean/@fams[r];
zpad[r]=Transpose[{zAlist[ℱ[r]][[21]],zDlist[ℱ[r]][[21]]}];
midp[r]=((zpad[r][[#[[1]]]]+zpad[r][[#[[2]]]])/2)&/@pars;,{r,10}];

Note that the pedigree is stored once, associated with 𝒩, not with the replicates 𝒩ℱ[r]. The correct calculation is to find the (Subscript[F, i,i]) which determines the variance released by meiosis within each diploid individual.  ft gives (Subscript[F, 12]) for each pair; (Subscript[F, 11]) is the probability of identity between random genes, (1+(Subscript[F, 11]))/2.  ft3 gives ((Subscript[F, 11])+(Subscript[F, 22]))/2.

Note that these quantities are defined only for 100 random pairs

ft=(ff=Fx[𝒩,20][{{{#[[1]],1},{#[[2]],2}}}];If[#[[1]]==#[[2]],((1+ff)/2),ff])&/@pars;
ft2=Table[Fx[𝒩,20][{{{i,1},{i,2}}}],{i,ninds}];
ft3=((ft2[[#[[1]]]]+ft2[[#[[2]]]])/2)&/@pars;

This was not stored & needs fresh evaluation

ft112=((Fx[𝒩,20][{{{#[[2]],1},{#[[2]],2},{#[[1]],1}}}]+Fx[𝒩,20][{{{#[[2]],1},{#[[2]],2},{#[[1]],1}}}])/2)&/@pars;

ft122=((Fx[𝒩,20][{{{#[[1]],1},{#[[1]],2},{#[[2]],1}}}]+Fx[𝒩,20][{{{#[[2]],1},{#[[2]],2},{#[[1]],1}}}])/2)&/@pars;
ft1122=Fx[𝒩,20][{{{#[[1]],1},{#[[1]],2},{#[[2]],1},{#[[2]],2}}}]&/@pars;
ft1212t=Fx[𝒩,20][{{{#[[1]],1},{#[[2]],1}},{{#[[1]],2},{#[[2]],2}}}]&/@pars;
ft1122t=Fx[𝒩,20][{{{#[[1]],1},{#[[1]],2}},{{#[[2]],1},{#[[2]],2}}}]&/@pars;

Save["Neutral identities 16th Jan",{ft,ft3,ft122,ft1122,ft1212t,ft1122t}];

((Subscript[σ, DI]^2)/4)(3(Subscript[F, 12])-(Subscript[F, 1122])-(Subscript[F, 112])-(Subscript[F, 122])) +((Subscript[V, D])/4)(3-(Subscript[F, 11])-(Subscript[F, 22])-(Subscript[F, 1122])-(Subscript[(OverTilde[F]), 1122])-(3/2)(Subscript[(OverTilde[F]), 1212]))

f1=(3ft-ft1122-2ft122)/4;f2=(3-2ft3-ft1122-ft1122t-(3/2)ft1212t)/4;
MinMax/@{σ2DI[pr,αh]f1,Vd[pr,αh]f2}

Working

fitTable=Table[Fit[Transpose[{f1,f2,vr[r][[All,2]]}],{1,F1,F2},{F1,F2}],{r,10}];
Join[fitTable,{Mean[fitTable]//Simplify,σ2DI[pr,αh]F1+Vd[pr,αh]F2}]//TableForm

GraphicsRow[{Histogram[fitTable/.{F1->0,F2->0},10],Histogram[D[fitTable,F1],10],Histogram[D[fitTable,F2],10]}]

vm=Mean[Table[vr[r][[All,2]],{r,1,10}]];

Show[ListPlot[Transpose[{σ2DI[pr,αh]f1+Vd[pr,αh]f2,vm}]],Plot[x,{x,0.05,0.07}],PlotRange->{{0,0.08},{0,0.08}},AxesOrigin->{0,0}]

Flatten[Table[Transpose[{σ2DI[pr,αh]f1,vr[r][[All,2]]}],{r,10}],1]//ListPlot

Final runs 9.2.2020

Strategic comments

The IM predicts a multivariate normal distribution of the individuals in the pedigree, with mean and covariance that depend on the various identity coefficients.  Because we run replicates across a fixed pedigree, sampling from a given reference population. these can be predicted from the theory.  One can present this, as below, by plotting the variances against (say) (Subscript[F, i,j]).  I expect that there will be a systematic error, because as the population becomes more inbred, all the coefficients increase, and so a regression on one of the coefficients will be steeper than expected from just that coefficient. The noise around the regression will be due to random variation in the other F, plus the inevitable sampling error arising from estimation of a variance from 1000 values.

A more rigorous approach would be to plot the simulated variance against the prediction, based on all the F’s.  The variance around this should be a (χ^2) distribution with ~1000 d.f.. One can then ask how much of the variation is accounted for by pairwise F, and how much better the prediction is when all the F are included.  This should be feasible, though calculating fourth order identities for all generations may be challenging: there are 27405 4 way combinations to store for each of the 50 generations, so I guess here will be ~ 2 million coefficients altogether.

I simulate 50 generations each with 30 individuals, and so one can test 1500 variances and 1124250 covariances.  That would be excessive, but one might at least check covariances within chosen generations.

One could also check the shape of the distribution - but I propose checking normality only for the families generated below.

Setting up

Restoring what’s needed for the neutral families

See above for how these data were generated

This gives the pedigrees and the initial population: {pedigree,pedigreeList,initialP,αh}

<<"Pedigrees N=30 neutral";<<"Initial population N=30 II";
pr=Mean[Flatten[initialP[[1]],1]]//N;

This sets up the families and the parent pairs: {fams,pars,zAlist,zDlist}

<<"Neutral families 16th Jan";(* 555Mb*)

It is quick to set up tables of family means and midparent values:

nreps =1000;ninds=30;nloci=1000;T=50;
Do[cvf[r]=Covariance/@fams[r];
mnf[r]=Mean/@fams[r];
zpad[r]=Transpose[{zAlist[ℱ[r]][[21]],zDlist[ℱ[r]][[21]]}];
midp[r]=((zpad[r][[#[[1]]]]+zpad[r][[#[[2]]]])/2)&/@pars;,{r,10}];

This sets up the identity coefficients: {ft,ft3,ft122,ft1122,ft1212t,ft1122t}. Note that the pedigree is stored once, associated with 𝒩, not with the replicates ℱ[r]. The correct calculation is to find the (Subscript[F, i,i]) which determines the variance released by meiosis within each diploid individual.  ft gives (Subscript[F, 12]) for each pair; (Subscript[F, 11]) is the probability of identity between random genes, (1+(Subscript[F, 11]))/2.  ft3 gives ((Subscript[F, 11])+(Subscript[F, 22]))/2.

<<"Neutral identities 16th Jan";

This calculates pairwise coefficients within and between individuals:

<<"Pairwise F N=30 neutral II";

pfw[t_]:=Table[Fx[𝒩,t][{{{i,1},{i,2}}}],{i,ninds}];
pfb[t_]:=Flatten[Table[Fx[𝒩,t][{{{i,1},{j,1}}}],{i,2,ninds},{j,1,i-1}]];

This recovers neutral simulation results {pList, zAlist, zDlist} for 𝒩[r], 𝒮[r]

<<"p, zA, zD, N=30 neutral III";

<<"p, zA, zD, N=30 selected k=5";

Initial values of (Subscript[V, a]),(Subscript[V, d]),cov[A,D],ℐ, (Subscript[σ, DI]^2), (Subscript[σ, ADI])

These are the initial components, measured in the reference population. Note that (SuperStar[ι])=(Subscript[V, D]) for biallelic loci.

{Va[pr,αh],Vd[pr,αh],Vad[pr,αh],ℐ[pr,αh],ℐ2[pr,αh],σADI[pr,αh],σ2DI[pr,αh] }

There is a subtle issue here, however: these will not be precisely the same as for genes sampled from the single actual base population. For example, 𝔼[(SuperStar[ζζ])] equals p q for identical genes, but equals -((p q)/(2N-1)) if two distinct genes are sampled from a finite population with known frequency p.  There will be some randomness due to having sampled a particular set of genes - though this will decline with the # of loci.

Predicting changes through time

Mean (OverBar[A]), (OverBar[D]) through time

This shows the mean (OverBar[A]),(OverBar[D]) through time, averaged over 1000 replicates. (OverBar[D]) is compared to ι (Subscript[(OverBar[F]), w]) (red).  The mean (OverBar[A]) fluctuates close to zero, as expected.  The mean (OverBar[D]) decreases as the population becomes more inbred.  There is a slight discrepancy from ℐ F which increases through time; the data fit (1-(1/N))F (green):

GraphicsRow[{ListLinePlot[Transpose[{Range[0,50],Mean[Table[Mean/@zAlist[𝒩[r]],{r,1000}]]}],AxesLabel->{"t","(OverBar[A])"}],Show[ListLinePlot[Transpose[{Range[0,50],Mean[Table[Mean/@zDlist[𝒩[r]],{r,1000}]]}]],
ListLinePlot[Table[{t,ℐ[pr,αh] Mean[pfw[t]]},{t,0,50}],PlotStyle->Red],ListLinePlot[Table[{t,(1-(1/30))ℐ[pr,αh] Mean[pfw[t]]},{t,0,50}],PlotStyle->Green],AxesLabel->{"t","(OverBar[D])"},PlotRange->All]}]

* Change in mean and variance (Fig 2)

Figure.  Changes of the mean and variance over 50 generations of neutral evolution.  The top row shows a single replicate, whilst the bottom row shows the average over 300 replicates.  The left column shows the means (G=A+D, A, D; black, blue, red), whilst the right column shows the variance components ((Subscript[V, G]), (Subscript[V, A]), (Subscript[V, D]), (Subscript[V, A,D]); black, blue, red, purple). On the right,  solid lines show the total variances and covariance, whilst the dashed lines show the genic component.  These differ through the contribution of linkage disequilibrium, which generates substantial variation.  The genic component changes smoothly, as expected with a large number (M=1000) of loci.   Simulations are made on a single pedigree with 30 individuals; variance components are measured relative to the base population.

nreps=300;T=50;c1={Black,Blue,Red};c2={Black,Blue,Red,Purple};
dash[c_]:={#,Dashing[{0.01,0.02}]}&/@c;
GraphicsGrid[{{fPlot[ℳ[1000,1,300],{zbGlist,zbAlist,zbDlist},c1,{1},T,{-1,1},"(OverBar[G]),(OverBar[A]),(OverBar[D])",Ticks->{{0,20,40},{-1,0,1}}],
Show[fPlot[ℳ[1000,1,300],{vGlist,vAlist,vDlist,vADlist},c2,{1},T,{-0.1,0.4},"(Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}],
fPlot[ℳ[1000,1,300],{vGGlist,vAGlist,vDGlist,vADGlist},dash[c2],{1},T,{-0.1,0.4},"(Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}]]},{fPlot[ℳ[1000,1,300],{zbGlist,zbAlist,zbDlist},c1,nreps,T,{-0.1,0.1},"(OverBar[G]),(OverBar[A]),(OverBar[D])",Ticks->{{0,20,40},{-0.1,0,0.1}}],
Show[fPlot[ℳ[1000,1,300],{vGlist,vAlist,vDlist,vADlist},c2,nreps,T,{-0.1,0.4},"(Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}],
fPlot[ℳ[1000,1,300],{vGGlist,vAGlist,vDGlist,vADGlist},dash[c2],nreps,T,{-0.1,0.4},"(Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}]]}}]

This is a version with just the bottom row, showing averages over 300 replicates

nreps=300;T=50;c1={Black,Blue,Red};c2={Black,Blue,Red,Purple};
dash[c_]:={#,Dashing[{0.01,0.02}]}&/@c;
GraphicsRow[{
fPlot[ℳ[1000,1,300],{zbGlist,zbAlist,zbDlist},c1,nreps,T,{-0.1,0.1},"(OverBar[G]),(OverBar[A]),(OverBar[D])",Ticks->{{0,20,40},{-0.1,0,0.1}}],
Show[fPlot[ℳ[1000,1,300],{vGlist,vAlist,vDlist,vADlist},dash[c2],nreps,T,{-0.1,0.4},"(Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}],
fPlot[ℳ[1000,1,300],{vGGlist,vAGlist,vDGlist,vADGlist},c2,nreps,T,{-0.1,0.4},"(Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}]]}]

* A, D vs F (Fig 3)

The relation between (Subscript[F, 12]) and A,D (left, right), in generations 5, 10, 20, 40 (black, blue, purple, red).  SInce the same pedigree is used for all replicates, each individual is characterised by a single (Subscript[F, 12]).  The plots show A, D, averaged over 1000 replicates, for each individual; thus, there are N=30 points in each generation.  There is no systematic effect of F on the mean A, whereas the mean D decreases as ℐ F = -0.53F (red line at right).

Note that the upward drift of A is not significant (see above): there is a strong correlation amongst genes within each replicate.

gr[f_,l_]:=Show[MapThread[ListPlot[Transpose[{pfw[#1],Mean[Table[f[𝒩[r]][[#1+1]],{r,nreps}]]}],PlotStyle->#2]&,{{5,10,20,40},{Black,Blue,Purple,Red}}],PlotRange->All,AxesLabel->{"(Subscript[F, 12])",l}];
GraphicsRow[{gr[zAlist,"A"],Show[gr[zDlist,"D"],Plot[ℐ[pr,αh] f,{f,0,0.55},PlotStyle->Red]]}]

Figure. The expected dominance deviation, 𝔼[D], declines with the identity between the two genes within each individual, (Subscript[F, 11]).  Each point is the mean D over 1000 replicates, made by dropping genes down a single pedigree. For generations 5, 10, 20, 40  (black, blue, purple, red), values are shown for each of the  N=30 individuals. The black line shows the prediction 𝔼[D]=ℐ (Subscript[F, 11]).

bs={FontFamily->"Times",FontSize->14};Show[MapThread[ListPlot[Transpose[{pfw[#1],Mean[Table[zDlist[𝒩[r]][[#1+1]],{r,nreps}]]}],PlotStyle->#2]&,{{5,10,20,40},{Black,Blue,Purple,Red}}],Plot[ℐ[pr,αh] f,{f,0,0.55},PlotStyle->Black],PlotRange->All,AxesLabel->{"(Subscript[F, 11])","𝔼[D]"},BaseStyle->bs,Ticks->{{0,0.2,0.4},{-0.3,-0.2,-0.1,0}}]

* Predicting var(A), var(D), cov(A,D) across individuals (Fig 4)

This constructs tables of the variance & covariance of A, D across replicates, for each of the 30 individuals, in a given generation:

varA[t_Integer]:=Variance[Table[zAlist[𝒩[r]][[t+1]],{r,nreps}]];
varD[t_Integer]:=Variance[Table[zDlist[𝒩[r]][[t+1]],{r,nreps}]];
varAD[t_Integer]:=Table[Covariance@@Transpose[Table[{zAlist[𝒩[r]][[t+1]],zDlist[𝒩[r]][[t+1]]},{r,nreps}]][[All,All,k]],{k,ninds}];

Figure. The variance components (Subscript[V, A]), (Subscript[V, D]), (Subscript[V, A,D]) are predicted by the variance between the two genes within each individual, (Subscript[F, 11]).  Variance components are shown for generations 5, 10, 20, 40 (black,… red); as in Fig. xx, there are N=30 points for each generation.

bs={FontFamily->"Times",FontSize->14};gr[f_,t_,c_]:=ListPlot[Transpose[{pfw[t],f[t]}],PlotStyle->c,BaseStyle->bs];
gr2[f_,l_,yr_,ty_,eq_]:=gr2[f,l,yr,ty,eq,""];
gr2[f_,l_,yr_,ty_,eq_,pl_]:=Show[MapThread[gr[f,#1,#2]&,{{5,10,20,40},{Black,Blue,Purple,Red}}],Plot[eq,{F,0,0.6},PlotStyle->Black],PlotRange->{{0,0.6},yr},AxesOrigin->{0,0},AxesLabel->{"(Subscript[F, 11])",l},BaseStyle->bs,Ticks->{{0,0.2,0.4,0.6},ty},PlotLabel->pl];
GraphicsColumn[{gr2[varA,"(Subscript[V, A])",{0,0.45},{0,0.2,0.4},Va[pr,αh](1+F),"(Subscript[V, A])=(Subscript[σ, A]^2)(1+(Subscript[F, 11]))"],
gr2[varD,"(Subscript[V, D])",{0,0.16},{0,0.05,0.1,0.15},(1-(F^2))Vd[pr,αh]+F σ2DI[pr,αh],"(Subscript[V, D])=(Subscript[σ, D]^2)(1-(Subscript[F, 11]^2))+(Subscript[σ, D I]^2)(Subscript[F, 11])"],
gr2[varAD,"(Subscript[V, A, D])",{-0.1,0.01},{-0.1,-0.05,0},F σADI[pr,αh],"(Subscript[V, A, D])=(Subscript[σ, A D I])(Subscript[F, 11])"]}]

As a row

bs={FontFamily->"Times",FontSize->14};gr[f_,t_,c_]:=ListPlot[Transpose[{pfw[t],f[t]}],PlotStyle->c,BaseStyle->bs];
gr2[f_,l_,yr_,ty_,eq_]:=gr2[f,l,yr,ty,eq,""];
gr2[f_,l_,yr_,ty_,eq_,pl_]:=Show[MapThread[gr[f,#1,#2]&,{{5,10,20,40},{Black,Blue,Purple,Red}}],Plot[eq,{F,0,0.6},PlotStyle->Black],PlotRange->{{0,0.6},yr},AxesOrigin->{0,0},AxesLabel->{"(Subscript[F, 11])",l},BaseStyle->bs,Ticks->{{0,0.2,0.4,0.6},ty},PlotLabel->pl];
GraphicsRow[{gr2[varA,"(Subscript[V, A])",{0,0.45},{0,0.2,0.4},Va[pr,αh](1+F),"(Subscript[V, A])=(Subscript[σ, A]^2)(1+(Subscript[F, 11]))"],
gr2[varD,"(Subscript[V, D])",{0,0.16},{0,0.05,0.1,0.15},(1-(F^2))Vd[pr,αh]+F σ2DI[pr,αh],"(Subscript[V, D])=(Subscript[σ, D]^2)(1-(Subscript[F, 11]^2))+(Subscript[σ, D I]^2)(Subscript[F, 11])"],
gr2[varAD,"(Subscript[V, A, D])",{-0.1,0.01},{-0.1,-0.05,0},F σADI[pr,αh],"(Subscript[V, A, D])=(Subscript[σ, A D I])(Subscript[F, 11])"]}]

As a row without the equations

bs={FontFamily->"Times",FontSize->14};gr[f_,t_,c_]:=ListPlot[Transpose[{pfw[t],f[t]}],PlotStyle->c,BaseStyle->bs];
gr2[f_,l_,yr_,ty_,eq_]:=gr2[f,l,yr,ty,eq,""];
gr2[f_,l_,yr_,ty_,eq_,pl_]:=Show[MapThread[gr[f,#1,#2]&,{{5,10,20,40},{Black,Blue,Purple,Red}}],Plot[eq,{F,0,0.6},PlotStyle->Black],PlotRange->{{0,0.6},yr},AxesOrigin->{0,0},AxesLabel->{"(Subscript[F, 11])",l},BaseStyle->bs,Ticks->{{0,0.2,0.4,0.6},ty},PlotLabel->pl];
GraphicsRow[{gr2[varA,"(Subscript[V, A])",{0,0.45},{0,0.2,0.4},Va[pr,αh](1+F)],
gr2[varD,"(Subscript[V, D])",{0,0.16},{0,0.05,0.1,0.15},(1-(F^2))Vd[pr,αh]+F σ2DI[pr,αh]],
gr2[varAD,"(Subscript[V, A, D])",{-0.1,0.01},{-0.1,-0.05,0},F σADI[pr,αh]]}]

Version without the equations

GraphicsColumn[{gr2[varA,"(Subscript[V, A])",{0,0.45},{0,0.2,0.4},Va[pr,αh](1+F)],
gr2[varD,"(Subscript[V, D])",{0,0.16},{0,0.05,0.1,0.15},(1-(F^2))Vd[pr,αh]+F σ2DI[pr,αh]],
gr2[varAD,"(Subscript[V, A, D])",{-0.1,0.01},{-0.1,-0.05,0},F σADI[pr,αh]]}]

Predictions for families

Prediction for residuals within families:

var[A]= ((Subscript[V, A])/2)(1-(((Subscript[F, 11])+(Subscript[F, 22]))/2)) 
cov[A,D]= ((Subscript[σ, ADI])/2)((Subscript[F, 12])-(((Subscript[F, 112])+(Subscript[F, 122]))/2)) 
var[D]= ((Subscript[σ, DI]^2)/4)(3(Subscript[F, 12])-(Subscript[F, 1122])-(Subscript[F, 112])-(Subscript[F, 122])) +((Subscript[V, D])/4)(3-(Subscript[F, 11])-(Subscript[F, 22])-(Subscript[F, 1122])-(Subscript[(OverTilde[F]), 1122])-(3/2)(Subscript[(OverTilde[F]), 1212]))

Prediction for variance of the shared component across families:

var[A]=((Subscript[V, A])/2)(1+(((Subscript[F, 11])+(Subscript[F, 22]))/2)+2(Subscript[F, 12])) +(Subscript[σ, ADI])((Subscript[F, 12])+(((Subscript[F, 112])+(Subscript[F, 122]))/2)) +((Subscript[σ, DI]^2)/4)((Subscript[F, 12])+(Subscript[F, 112])+(Subscript[F, 122])+(Subscript[F, 1122])) +((Subscript[V, D])/4) (1 + (Subscript[F, 11])+(Subscript[F, 22])+(Subscript[F, 1122])+(Subscript[(OverTilde[F]), 1122])+(3/2)(Subscript[(OverTilde[F]), 1212])-4(Subscript[F, 12]^2))

(Subscript[V, w]) vs ((Subscript[F, ii])+(Subscript[F, jj]))/2

These tables are for plots of (Subscript[V, W]) against (1+((Subscript[F, ii])+(Subscript[F, jj]))/2)

res[r_]:=Transpose/@Table[fams[r][[pr,All,k]]-mnf[r][[pr,k]],{pr,1,100},{k,1,2}];
vr[r_]:=Variance/@res[r];
fvt[r_,k_]:=Transpose[{ft3,vr[r][[All,k]]}];
fit[r_,k_]:=Fit[fvt[r,k],{1,f},f];
ftt[k_]:=fit[#,k]&/@Range[10];

This shows 10 replicate fits, each based on 100 families of 1000 offspring.

TableForm[Join[ftt[#],{Simplify[Mean[ftt[#]]],Va[pr,αh](1-f)/2}]]&/@{1,2}

These tables are for plots of (Subscript[V, W]) against the midparent trait value:

fvtA[r_,k_]:=Transpose[{midp[r][[All,k]],vr[r][[All,k]]}];
fitA[r_,k_]:=Fit[fvtA[r,k],{1,f},f];
fttA[k_]:=fitA[#,k]&/@Range[10];

The variance of A within families does not depend on the mid-parental value of A (top left), but decreases with (Subscript[V, A])(1-((Subscript[F, 11])+(Subscript[F, 22]))/2) (red line at top right); the blue line is the least-squares regression. The second row shows the same for the variance of D within families; this does not depend on the mid-parental D (bottom left). Need to work out the prediction for (Subscript[V, D]) and worry about whether I should fit the more complicated formula with all the F’s  These plots superimpose ten replicates, each with 100 families of 1000 offspring.

mnftt=Mean[ftt[1]];mnfttA=Mean[fttA[1]];mnftt2=Mean[ftt[2]];mnfttA2=Mean[fttA[2]];GraphicsGrid[{{Show[ListPlot[fvtA[#,1]]&/@Range[10],AxesLabel->{"((Subscript[(OverBar[A]), i])+(Subscript[(OverBar[A]), j]))/2","(Subscript[V, A, WF])"},PlotRange->{{-1.5,1.5},{0,0.15}},AxesOrigin->{0,0}],
Show[ListPlot[fvt[#,1]]&/@Range[10],Plot[{mnftt,Va[pr,αh](1-f)/2},{f,0.2,0.4},PlotStyle->{Blue,Red}],AxesLabel->{"((Subscript[F, ii])+(Subscript[F, jj]))/2","(Subscript[V, A, WF])"},PlotRange->{{0.2,0.35},{0,0.15}},AxesOrigin->{0.2,0}]},
{Show[ListPlot[fvtA[#,2]]&/@Range[10],AxesLabel->{"((Subscript[(OverBar[D]), i])+(Subscript[(OverBar[D]), j]))/2","(Subscript[V, D, WF])"},PlotRange->{{-1.5,1.5},{0,0.08}},AxesOrigin->{0,0}],
Show[ListPlot[fvt[#,2]]&/@Range[10],Plot[{mnftt2,Va[pr,αh](1-f)/2},{f,0.2,0.4},PlotStyle->{Blue,Red}],AxesLabel->{"((Subscript[F, ii])+(Subscript[F, jj]))/2","(Subscript[V, D, WF])"},PlotRange->{{0.2,0.35},{0,0.08}},AxesOrigin->{0.2,0}]}}]
{mnfttA,mnfttA2}

** (Co)variance of A,D within families (Fig 5)

Clear[cvt]

These tables are for plots of (Subscript[V, W]) against (1+((Subscript[F, ii])+(Subscript[F, jj]))/2)

res[r_]:=Transpose/@Table[fams[r][[pr,All,k]]-mnf[r][[pr,k]],{pr,1,100},{k,1,2}];
vr[r_]:=Variance/@res[r];
cv[r_]:=Covariance@@Transpose[#]&/@res[r];
fvt[r_,k_]:=Transpose[{ft3,vr[r][[All,k]]}];
cvt[r_]:=Transpose[{(ft-ft122),cv[r]}];(* prediction for cov[A,D] *)
vdt1[r_]:=Transpose[{(3ft-ft1122-ft112-ft122),vr[r][[All,2]]}];(* prediction for var[D] *)
vdt2[r_]:=Transpose[{(3-2ft3-ft1122-ft1122t-(3/2)ft1212t),vr[r][[All,2]]}];(* prediction for var[D] *)
fit[r_,k_]:=Fit[fvt[r,k],{1,f},f];
ftt[k_]:=fit[#,k]&/@Range[10];

ft112

a=(3ft-ft1122-ft112-ft122);
b=(3-2ft3-ft1122-ft1122t-(3/2)ft1212t);
y=Mean[vr[#][[All,2]]&/@Range[10]]

tt=Transpose[{a,b,y}];tts=Pick[tt,(#≤0.5)&/@ft];
vdf=Fit[tts,{x,z},{x,z}];
{vdf,σ2DI[pr,αh]/4,Vd[pr,αh]/4}

Show[ListPlot[Transpose[{a,y}]],Plot[(1/4)σ2DI[pr,αh]f+(1/4)Vd[pr,αh]Mean[b],{f,0.4,0.8},PlotStyle->Red],PlotRange->{{0,0.9},{0,0.08}},AxesOrigin->{0,0}]

Show[ListPlot[Transpose[{(1/4)(σ2DI[pr,αh]a+Vd[pr,αh]b),y}]],Plot[x,{x,0.05,0.1},PlotStyle->Red],PlotRange->{{0,0.08},{0,0.08}},AxesOrigin->{0,0},AxesLabel->{"(Subscript[σ, D]^2), pred.","(Subscript[σ, D]^2)"}]

xx=(1/4)(σ2DI[pr,αh]a+Vd[pr,αh]b);
tt=Transpose[{xx,y}];tts=Pick[tt,(#≤0.5)&/@ft];
vdf=Fit[tts,{x},{x}];
Show[ListPlot[Transpose[{xx,vr[#][[All,2]]}]]&/@Range[10],Plot[{x,vdf},{x,0.05,0.1},PlotStyle->{Red,Blue}],PlotRange->{{0,0.08},{0,0.08}},AxesOrigin->{0,0},AxesLabel->{"(Subscript[σ, D]^2), pred.","(Subscript[σ, D]^2)"}]

Show[ListPlot[Transpose[{b,y}]],Plot[(1/4)σ2DI[pr,αh]Mean[a]+(1/4)Vd[pr,αh]f,{f,1,2.5},PlotStyle->Red],PlotRange->All]

GraphicsRow[Histogram[#,20]&/@{a,b}]

var[A]= ((Subscript[V, A])/2)(1-(((Subscript[F, 11])+(Subscript[F, 22]))/2)) 
cov[A,D]= ((Subscript[σ, ADI])/2)((Subscript[F, 12])-(((Subscript[F, 112])+(Subscript[F, 122]))/2)) 
var[D]= ((Subscript[σ, DI]^2)/4)(3(Subscript[F, 12])-(Subscript[F, 1122])-(Subscript[F, 112])-(Subscript[F, 122])) +((Subscript[V, D])/4)(3-(Subscript[F, 11])-(Subscript[F, 22])-(Subscript[F, 1122])-(Subscript[(OverTilde[F]), 1122])-(3/2)(Subscript[(OverTilde[F]), 1212]))

mftt1=Mean[ftt[1]];mftt2=Mean[ftt[2]];
xx=(1/4)(σ2DI[pr,αh]a+Vd[pr,αh]b);
tt=Transpose[{xx,y}];tts=Pick[tt,(#≤0.5)&/@ft];
vdf=Fit[tts,{x},{x}];

Figure.  The variances and covariance within families between the residual additive and dominance deviations, R, S.  The blue lines shows a least-square regression, and the red lines, the prediction (Eq. xx).  100 pairs of parents were chosen at random, and from each, 1000 offspring were generated. This was repeated for 10 replicates, made using the same pedigree and the same set of parents; within-family variances were averaged over replicates.  In each figure, each 100 point corresponds to one of the 100 pairs of parents. The five outlier points at right are families produced by selfing.

The predictions are (Subscript[V, R])=(Subscript[σ, A]^2)(1-(Subscript[F, W]))/2, where (Subscript[F, W])=((Subscript[F, 11])+(Subscript[F, 22]))/2, the within-individual identity averaged over parents 1,2;  (Subscript[V, S]) =(Subscript[σ, D I]^2)(3(Subscript[F, 12])-(Subscript[F, 1122])-(Subscript[F, 112])-(Subscript[F, 122]))/4+(Subscript[σ, D]^2)(3-(Subscript[F, 11])-(Subscript[F, 22])-(Subscript[F, 1122])-(Subscript[(OverTilde[F]), 1122])-(3/2)(Subscript[(OverTilde[F]), 1212]))/4; (Subscript[V, R,S])=(Subscript[σ, A D I])((Subscript[F, 12])-(Subscript[F, (3)]))/2, where (Subscript[F, (3)])=((Subscript[F, 112])+(Subscript[F, 122]))/2.

bs={FontFamily->"Times",FontSize->14};GraphicsRow[{Show[ListPlot[Mean[fvt[#,1]&/@Range[10]]],Plot[{mftt1,Va[pr,αh](1-f)/2},{f,0.2,0.35},PlotStyle->{Blue,Red}],AxesLabel->{"1-(Subscript[F, w])","(Subscript[V, R])"},PlotRange->{{0.2,0.35},{0,0.12}},Ticks->{{0.2,0.25,0.3,0.35},{0,0.05,0.1}},AxesOrigin->{0.2,0},BaseStyle->bs],(*Show[ListPlot[fvt[#,2]]&/@Range[10],Plot[{mftt2,Va[pr,αh](1-f)/2},{f,0.2,0.4},PlotStyle->{Blue,Red}],AxesLabel->{"((Subscript[F, 11])+(Subscript[F, 22]))/2","(Subscript[σ, D, WF]^2)"},PlotRange->{{0.2,0.35},{0,0.08}},AxesOrigin->{0.2,0}]*)
Show[ListPlot[Mean[Transpose[{xx,vr[#][[All,2]]}]&/@Range[10]]],Plot[{x,vdf},{x,0.05,0.1},PlotStyle->{Red,Blue}],PlotRange->{{0.04,0.08},{0.04,0.08}},Ticks->{{0.04,0.06,0.08},{0.04,0.06,0.08}},AspectRatio->1,AxesOrigin->{0.04,0.04},AxesLabel->{"(Subscript[V, S]), pred.","(Subscript[V, S])"},BaseStyle->bs],
Show[ListPlot[Mean[cvt/@Range[10]]],Plot[{mftt2,σADI[pr,αh]df/2},{df,0.1,0.4},PlotStyle->{Blue,Red}],AxesLabel->{"(Subscript[F, 12])-(Subscript[F, (3)])","(Subscript[V, R, S])"},Ticks->{{0.2,0.3,0.4},{-0.03,-0.02,-0.01,0}},PlotRange->{{0.1,0.4},{-0.03,0}},AxesOrigin->{0.1,0},BaseStyle->bs]}]

Version showing each of the 10 replicates, rather than the mean: 1000 points rather than 100

bs={FontFamily->"Times",FontSize->14};GraphicsRow[{Show[ListPlot[fvt[#,1]]&/@Range[10],Plot[{mftt1,Va[pr,αh](1-f)/2},{f,0.2,0.4},PlotStyle->{Blue,Red}],AxesLabel->{"1-(Subscript[F, w])","(Subscript[V, R])"},PlotRange->{{0.2,0.35},{0,0.15}},Ticks->{{0.2,0.25,0.3},{0,0.05,0.1}},AxesOrigin->{0.2,0},BaseStyle->bs],(*Show[ListPlot[fvt[#,2]]&/@Range[10],Plot[{mftt2,Va[pr,αh](1-f)/2},{f,0.2,0.4},PlotStyle->{Blue,Red}],AxesLabel->{"((Subscript[F, 11])+(Subscript[F, 22]))/2","(Subscript[σ, D, WF]^2)"},PlotRange->{{0.2,0.35},{0,0.08}},AxesOrigin->{0.2,0}]*)
Show[ListPlot[Transpose[{xx,vr[#][[All,2]]}]]&/@Range[10],Plot[{x,vdf},{x,0.05,0.1},PlotStyle->{Red,Blue}],PlotRange->{{0,0.08},{0,0.08}},Ticks->{{0,0.04,0.08},{0,0.04,0.08}},AspectRatio->1,AxesOrigin->{0,0},AxesLabel->{"(Subscript[V, S]), pred.","(Subscript[V, S])"},BaseStyle->bs],
Show[ListPlot[cvt[#]]&/@Range[10],Plot[{mftt2,σADI[pr,αh]df/2},{df,0.1,0.4},PlotStyle->{Blue,Red}],AxesLabel->{"(Subscript[F, 12])-(Subscript[F, (3)])","(Subscript[V, R, S])"},Ticks->{{0.2,0.3,0.4},{-0.02,-0.01,0,0.01}},PlotRange->{{0.1,0.4},{-0.03,0.01}},AxesOrigin->{0.1,0},BaseStyle->bs]}]

Selection

Setting up replicates with different # of loci

Read in the pedigree, 𝒩:

<<"Pedigrees N=30 neutral";

This sets up a WF population (but with no selfing) of N=30 diploid individuals, and varying # of loci, M. Dominance is h=0,1; probability of h=1 is (1/2)+(2/(Sqrt[M])), which should ensure (Subscript[V, D]) of order 1.  Allele frequencies are sampled from a large base population with allele frequencies distributed as (1/(p q)).  The reference frequencies are those in the actual sampled population.

I set up a single array of α, h and p for M=(10^4), and then choose subsets of this, scaled so that the total genetic variance in the base population is the same, at 0.34 (as in the examples above). However, (Subscript[V, D]) may vary somewhat by chance when M is small.  I will use M=100, 300, 1000, 3000, (10^4), and select for 50 generations.

ninds=30;nlMax=10000;nreps=300;

αhMax=RandomSample[setαh[nlMax,2.]]//N;
jr=randomJ[2ninds,0.5,0.999,nlMax];prMax=jr/(2ninds);
prM[M_]:=prM[M]=prMax[[1;;M]];
αhM[M_]:=αhM[M]=Module[{vm=Vg[prMax,αhMax],vv=Vg[prMax[[1;;M]],αhMax[[1;;M]]]},
({(Sqrt[(vm/vv)]),1}#)&/@αhMax[[1;;M]]];

Timing[(mm=#;makeSel[{𝒩,ℳ},mm,2,300,30,50];)&/@{100,300,1000,3000};]

Timing[makeSel[{𝒩,ℳ},10000,1,300,30,50];]

Timing[makeSel[{𝒩,ℳ},10000,2,300,30,50];]

Timing[(mm=#;makeSel[{𝒩,ℳ},mm,1,300,30,50];)&/@{100,300,1000,3000};]

Save["selected reps N=30 selected k=1,2 300 reps III", {prM,αhM,zAlist,zDlist,zbAlist,zbDlist,vAlist,vDlist,vADlist,vAGlist,vDGlist,vADGlist}];

<<"selected reps N=30 selected k=1,2 300 reps III";

Initial values of (Subscript[V, a]),(Subscript[V, d]),cov[A,D],ℐ,(SuperStar[ℐ]), (Subscript[σ, DI]^2), (Subscript[σ, ADI])

These are the initial components, measured in the reference population, for the various M. Note that (SuperStar[ι])=(Subscript[V, D]) for biallelic loci. The initial α were scaled to ensure the same total genic variance:

Prepend[{#,Va[prM[#],αhM[#]]+Vd[prM[#],αhM[#]],Va[prM[#],αhM[#]],Vd[prM[#],αhM[#]],Vad[prM[#],αhM[#]],ℐ[prM[#],αhM[#]],ℐ2[prM[#],αhM[#]],σADI[prM[#],αhM[#]],σ2DI[prM[#],αhM[#]] }&/@{100,300,1000,3000},{"M","(Subscript[V, g])","(Subscript[V, a])","(Subscript[V, d])","cov[A,D]","ℐ","(SuperStar[ℐ])"," (Subscript[σ, ADI])"," (Subscript[σ, DI]^2)"}]//TableForm

There is a subtle issue here, however: these will not be precisely the same as for genes sampled from the single actual base population. For example, 𝔼[(SuperStar[ζζ])] equals p q for identical genes, but equals -((p q)/(2N-1)) if two distinct genes are sampled from a finite population with known frequency p.  There will be some randomness due to having sampled a particular set of genes - though this will decline with the # of loci.

* Change in mean and variance (Fig 6, 7)

All components are measured relative to the base population. Top row: means across 1000 replicates; bottom row: 10 replicates, plotted individually. Left: mean (black), additive component (grey), dominance component (blue); right: actual variance (including LD); (Subscript[V, A]) (grey); (Subscript[V, D]) (blue); cov(A,D) (purple); (Subscript[V, A])+(Subscript[V, D])+2cov[A,D] (red).   Examining individual replicates (bottom row) we see that indeed, the actual (Subscript[V, G])(including LD) is much more variable than the underlying genic components.  I should reconstruct the genic components of variance (which would be much smoother). But, that requires knowing the allele frequencies.

Selection this strong causes an immediate drop in variance. However, this may just be due to LD. I need to calculate the genic variance, and also do runs with different numbers of loci.

Note: there were only 300 neutral replicates, but 1000 selected replicates

Figure. Comparisons between neutral vs. selected simulations (dashed, solid lines).  Top: change in means relative to initial value ( (OverBar[G])=(OverBar[A+D]), (OverBar[A]), (OverBar[D]); black, blue, red); middle: variances, including LD (var(A+D), var(A), var(D), cov(A,D); black, blue, red, purple).  Values are averages over 300 replicates, made with the same pedigree. There are  M=1000 loci. Selection is made within families: for each pair of parents, two individuals are generated, and the one with the largest value is selected.

nreps=300;T=50;c1={Black,Blue,Red};c2={Black,Blue,Red,Purple};
dash[c_]:={#,Dashing[{0.01,0.02}]}&/@c;
GraphicsColumn[{
Show[fPlot[ℳ[1000,#,300],{zbGlist,zbAlist,zbDlist},If[#==2,c1,dash[c1]],nreps,T,{-1.5,6},"(OverBar[G]),(OverBar[A]),(OverBar[D])",Ticks->{{0,20,40},{0,2,4,6}}]&/@{1,2}],
Show[fPlot[ℳ[1000,#,300],{vGlist,vAlist,vDlist,vADlist},If[#==2,c2,dash[c2]],nreps,T,{-0.1,0.4},"(Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}]&/@{1,2}],
Show[fPlot[ℳ[1000,#,300],{vGGlist,vAGlist,vDGlist,vADGlist},If[#==2,c2,dash[c2]],nreps,T,{-0.1,0.4},"genic (Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}]&/@{1,2}]}]

As a row

nreps=300;T=50;c1={Black,Blue,Red};c2={Black,Blue,Red,Purple};
dash[c_]:={#,Dashing[{0.01,0.02}]}&/@c;
GraphicsRow[{
Show[fPlot[ℳ[1000,#,300],{zbGlist,zbAlist,zbDlist},If[#==2,c1,dash[c1]],nreps,T,{-1.5,6},"(OverBar[G]),(OverBar[A]),(OverBar[D])",Ticks->{{0,20,40},{0,2,4,6}}]&/@{1,2}],
Show[fPlot[ℳ[1000,#,300],{vGlist,vAlist,vDlist,vADlist},If[#==2,c2,dash[c2]],nreps,T,{-0.1,0.4},"(Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}]&/@{1,2}],
Show[fPlot[ℳ[1000,#,300],{vGGlist,vAGlist,vDGlist,vADGlist},If[#==2,c2,dash[c2]],nreps,T,{-0.1,0.4},"genic (Subscript[V, G]),(Subscript[V, A]),(Subscript[V, D]),(Subscript[V, A, D])",Ticks->{{0,20,40},{0,0.2,0.4}}]&/@{1,2}]}]

Figure.   Convergence of the variance components at 50 generations to the infinitesimal limit, as the number of loci increases from M=100 to (10^4).   Simulations with 50% truncation selection are compared  with neutral simulations (solid, dashed lines).  Replicate simulations were generated as in Fig. xx, above.

locs={100,300,1000,3000,(10^4)};dash[c_]:={#,Dashing[{0.01,0.02}]}&/@c;
bs={FontFamily->"Times",FontSize->14};
tbl[f_,k_]:={#,Last[Mean[Table[f[ℳ[#,k,nreps][r]],{r,nreps}]]]}&/@locs;
gr[f_,c_,yr_,yl_]:=Show[ListLogLinearPlot[tbl[f,1],Joined->True,PlotStyle->dash[{c}],BaseStyle->bs,PlotRange->{{100,(10^4)},yr},AxesOrigin->{100,0},Ticks->{{100,1000,{(10^4),"(10^4)"}},{-0.02,0,0.02,0.04,0.06}}],
ListLogLinearPlot[tbl[f,2],Joined->True,PlotStyle->c],AxesLabel->{"M",yl}];GraphicsGrid[{{gr[vGGlist,Black,{0,0.07},"(Subscript[V, G])"],gr[vAGlist,Blue,{0,0.07},"(Subscript[V, A])"]},{gr[vDGlist,Red,{0,0.07},"(Subscript[V, D])"],gr[vADGlist,Purple,{-0.025,0},"(Subscript[V, A, D])"]}}]

Figure.   Convergence of the variance components at 50 generations to the infinitesimal limit, as the number of loci increases from M=100 to (10^4).   Points show the absolute difference between selected and neutral simulations, whilst lines show the least square fit.  The slopes are -0.62, -0.72, -0.70, -0.66 for (Subscript[V, G]), (Subscript[V, A]), (Subscript[V, D]), (Subscript[V, A,D]). Because the deviation of the covariance (Subscript[V, A,D]) changes sign (Fig. xx, bottom right), the fit for that case is for M=1000, 3000, (10^4). These slopes indicate that convergence is somewhat faster than with (M^(-1/2)).

tbl[f_]:={#,Last[Mean[Table[Abs[f[ℳ[#,2,nreps][r]]-f[ℳ[#,1,nreps][r]]],{r,nreps}]]]}&/@locs;
fitM[f_]:=fitM[f,0];fitM[f_,tk_]:=Fit[Log[Drop[tbl[f],tk]],{1,m},m];gr2[f_,c_,yl_]:=gr2[f,c,yl,0];gr2[f_,c_,yl_,tk_]:=Module[{ff=fitM[f,tk]},Show[ListLogLogPlot[tbl[f],PlotStyle->{c,PointSize[0.03]},BaseStyle->bs,PlotRange->{{90,1.1 (10^4)},{0.001,0.1}},Ticks->{{100,1000,{(10^4),"(10^4)"}},{0.001,0.01,0.1}},AxesOrigin->{100,0.001},AxesLabel->{"M",yl}],LogLogPlot[Exp[ff]/.m->Log[x],{x,100,(10^4)}]]];
GraphicsGrid[{
{gr2[vGGlist,Black,"(Subscript[V, G])"],gr2[vAGlist,Blue,"(Subscript[V, A])"]},{gr2[vDGlist,Red,"(Subscript[V, D])"],gr2[vADGlist,Purple,"(Subscript[V, A, D])",2]}}]
MapThread[fitM,{{vGGlist,vAGlist,vDGlist,vADGlist},{0,0,0,2}}]

Definitions

Identity coefficients

pedigreeList

pedigreeList::usage="pedigreeList[𝒫] stores a list of pedigrees associated with the symbol 𝒫.  This allows functions such as Fx to refer to the list without storing it.  It needs to be defined before use.  ** Here, a pedigree is represented as {{(Subscript[par, 1]),(Subscript[par, 2])},…}, using Parents[]";

genesToParents

genesToParents::usage="genesToParents[{{{1,1},{3,2}},…},{{3,4},…}] gives sets of IBD genes that point to the diploid parents. The first argument lists sets of genes that are IBD within, but not between, sets. The second argument lists parents, numbered consecutively. The result is expressed as terms of the form {{4,5},…}, which are weighted equally. ";

genesToParents[g:{___List},pars:{{_Integer,_Integer}..}]:=
If[OrderedQ[g]∧(And@@(OrderedQ/@g)),
Module[{inds,inds1,inds2,h,f,n1,n2},
inds=First/@Sort[Flatten[g,1]];
inds1=Select[Union[inds],Count[inds,#]==1&];
inds2=Select[Union[inds],Count[inds,#]==2&];
n1=Length[inds1];n2=Length[inds2];
h=Flatten[Outer[List,
Sequence@@pars[[inds1]],Sequence@@ConstantArray[{1,2},n2]],n1+n2-1];
f[s:{___Integer}]:=Join[
MapThread[{#1,_}->#2&,{inds1,Take[s,n1]}],
Flatten[MapThread[{{#1,1}->pars[[#1,#2]],{#1,2}->pars[[#1,3-#2]]}&,{inds2,Take[s,-n2]}],1]];Replace[g,f[#],{2}]&/@h],
genesToParents[Sort[Sort/@g],pars](* put into ordered form *)];

parentsToGenes

parentsToGenes::usage="parentsToGenes[{{1,3,5},{2,5}}] gives the joint probability of identity between genes drawn randomly from the specified individuals, as a list {{{…},w},…}.  parentsToGenes[{{{1,3,5},{2,5}},w}] returns a weighted list. For each distinct individual, genes are partitioned randomly into two sets, and the corresponding sets are merged if genes come from the same set. ";

parentsToGenes[g:{___List}]:=Module[{inds,sets,spl},
inds=Union[Flatten[g]];
sets=(First/@Position[g,#])&/@inds;
spl=makeAllSplits[sets];
{Replace[gatherSets[#[[1]]],{j_Integer,k:(1|2)}:>{inds[[j]],k},{2}],#[[2]]}&/@spl];

distinctQ

distinctQ::usage="distinctQ[{{a,b},{c,d}…}] tests whether the elements {a,b,c,d …} are distinct";

distinctQ[s:{___List}]:=Module[{sl},
sl=Flatten[s,1];Length[sl]==Length[Union[sl]]];

Fx

Fx::usage="Fx[𝒫,t][{{{1,1},{3,2}},…}] stores the probability that the listed sets of genes are all IBD within sets, but distinct from each other. This assumes that pedigreeList[𝒫] lists of parents for each individual, numbered consecutively, for each generation. The result is expressed in terms of a sum over Fx[𝒫,t][{{{4,1},{5,2}},{…},…},…}], where the sets point to genes in diploid individuals. Fx[𝓅,0][] is set to 0 or 1 as appropriate";

Fx[_,0][g:{___List}]:=
If[distinctQ[g],
If[Or@@(Length[#]>1&/@g),0.,1.],
Message[Fx::multipleGenes,g];$Failed];
(* If the sets are all distinct individuals, and all of length 1 (or 0), then F is 1. The initial F is set to a numerical value to avoid Mathematica getting bogged down in exact arithmetic *)

Fx[𝒫_,t_Integer][g:{___List}]:=Fx[𝒫,t][g]=Module[{g1},
If[distinctQ[g],
If[g=={},1,
If[OrderedQ[g]∧(And@@(OrderedQ/@g)),
Fx[𝒫,t][g]=
iterF[Fx[𝒫,t-1],g,pedigreeList[𝒫][[t]]],
Fx[𝒫,t][Sort[Sort/@g]](* put into ordered form *)]],
(* if genes appear more than once within sets, replace by the union *)
g1=Union/@g;
(* if the same gene appears in multiple sets, then F=0 *)
If[distinctQ[g1],Fx[𝒫,t][g1],0]
]];

Fx::multipleGenes="Genes appear more than once in `1`";

iterF

iterF::usage="iterF[f,{{{1,1},{3,2}},…},{{1,3},…}] iterates IBD amongst sets of genes, given a list of parents of diploid individuals, returning the result as an expression in the form ∑ w*f[{{{1,1},…},…}]. Fx sets f to Fx[𝒫,t-1]";

iterF[f_,g:{___List},pars:{{_Integer,_Integer}..}]:=Module[{gtp,tgtp,ptg},
gtp=genesToParents[g,pars];
(* The original version is inefficient, because gtp contains many equivalent sets.
Total[(((#[[2]])/(Length[gtp]))f[#[[1]]])&/@Flatten[parentsToGenes/@gtp,1]*)
tgtp=Tally[Map[Sort,gtp,2]];(* tally equivalent permutations *)
Total[(((#[[2]])/(Length[gtp]))f[#[[1]]])&/@Flatten[(ptg=parentsToGenes[#[[1]]];Transpose[{ptg[[All,1]],ptg[[All ,2]]#[[2]]}])&/@tgtp,1]]
(* for each permutation generated by genesToParents, produce the list of IBD sets, and weight by the multiplicity of that permutation *) ];

Understanding the algorithm

This sets up lists of the parents, over 50 generations, for 10 individuals:

pedigreeList[𝒫]=Table[Parents[randomPedigree[{10,10}]],{50}];

TableForm[pedigreeList[𝒫][[1;;3]],TableDepth->2]

Table[{Fx[𝒫,5][{{{k,1},{k,2}}}],Fx[𝒫,5][{{{k,1},{k,2}},{{j,1}}}]+Fx[𝒫,5][{{{k,1},{k,2},{j,1}}}]},{j,2,10},{k,1,j-1}]

Fx[𝒫,3]/@{{{{3,1},{3,2},{2,1},{2,2}}},{{{3,1},{3,2},{2,1},{2,2}}},{{{3,1},{2,1}},{{3,2},{2,2}}}}

pars=pedigreeList[𝒫][[3]];g={{{3,1},{3,2},{2,1},{2,2}}};
gtp=genesToParents[g,pars];
tgtp=Tally[Map[Sort,gtp,2]]

xx=(ptg=parentsToGenes[#[[1]]];Transpose[{ptg[[All,1]],ptg[[All ,2]]#[[2]]}])&/@tgtp

Total[(((#[[2]])/(Length[gtp]))f[#[[1]]])&/@Flatten[xx,1]]

iterF[f,g,pars]

pars=pedigreeList[𝒫][[3]];g1={{{3,1},{3,2}},{{2,1},{2,2}}};
gtp1=genesToParents[g1,pars];
tgtp1=Tally[Map[Sort,gtp1,2]]

Total[(((#[[2]])/(Length[gtp1]))f[#[[1]]])&/@Flatten[(ptg=parentsToGenes[#[[1]]];Transpose[{ptg[[All,1]],ptg[[All ,2]]#[[2]]}])&/@tgtp1,1]]

makeAllSplits

makeAllSplits::usage="makeAllSplits[{{1,2},…,{3,4}}] takes a list of the sets that descend from each individual, and lists all the ways that these can be validly assigned, together with the weights, (Subscript[w, i]).  The result is a list in the form {{{{{1},{2}},{{},{3}},…},(Subscript[w, 1])},…}, where the first entry lists the sets which descend from the two genes in each individual. ";

makeAllSplits[sets_List]:=Module[{aa,bb,ts,wts},
aa=(ts=Tally[Sort/@findSplit[#]];ts/.{a_List,j_Integer}:>{Union/@a,j (2^(-Length[#]))})&/@sets;
bb=Flatten[Outer[List,Sequence@@aa,1],Length[aa]-1];
wts=(Times@@#[[All,2]])&/@bb;
Transpose[{bb[[All,All,1]],wts}]];

For example, suppose that we have four sets of IBD genes, two of which descend from one individual, and two from the other.  This can only happen if 1 & 2 descend from the different genes in one individual, and 3 & 4 from the different genes in the other, which has probability 1/4:

makeAllSplits[{{1,2},{3,4}}]

gatherSets

gatherSets::usage="gatherSets[{{{1},{2}},…,{{},{3}}}] takes a list of the sets of genes that descend from each individual (Subscript[i, 1]),…,(Subscript[i, n]), and gathers together sets of genes that are IBD. An entry {{x},{y}} for individual a corresponds to genes {a,1} in set x and {a,2} in set y, whilst {{},{x}} corresponds to a gene {a,1} in set x.";

gatherSets[g:{{{___Integer},{_Integer}}...}]:=Module[{n=Length[g],xx,ss},
xx=Flatten[Table[g[[j]]/.{{{},{s_Integer}}:>{{s,{j,1}}},{{s1_Integer},{s2_Integer}}:>{{s1,{j,1}},{s2,{j,2}}}},{j,n}],1];
ss=Union[First/@xx];(Cases[xx,{#,_}]&/@ss)[[All,All,2]]];

findSplit

findSplit::usage="findSplit[{1,2,2}] finds valid partitions into two sets, out of the total (2^n)";

findSplit[x:{___Integer}]:=Select[{#,TakeAway[x,#]}&/@Subsets[x],(Length[Union[#[[1]]]]≤1)∧(Length[Union[#[[2]]]]≤1)&];

randomPedigree

randomPedigree::usage="randomPedigree[{n,p}] makes a random WF pedigree (but with no selfing), with n offspring and p parents. randomPedigree[{{(Subscript[n, 1]),(Subscript[n, 2]),…,(Subscript[n, k])},{(Subscript[p, 1]),…,(Subscript[p, k])}},{{(Subscript[m, 1, 1]),(Subscript[m, 1, 2]),…},…}] makes a random pedigree for the k-island model, with (Subscript[n, i]) offspring, (Subscript[p, i]) parents, and a fraction (Subscript[m, i, j]) in deme i coming from deme j. randomPedigree[{n,p},{{(Subscript[m, 1, 1]),(Subscript[m, 1, 2]),…},…}] assumes that all k demes have n offspring and p parents. There must be at least two parents per deme to avoid selfing. Returns a SparseArray. Selfing->True allows selfing, giving the WF model - only defined for one deme as yet";

Selfing::usage="Selfing is an option for randomPedigree that allows selfing";

Options[randomPedigree]={Selfing->False};

randomPedigree[{n_Integer,p_Integer},opts___Rule]:=Module[{c=ConstantArray[0,p],rp},
If[Selfing/.{opts}/.Options[randomPedigree],
SparseArray[Table[rp=RandomInteger[{1,p},{2}];
If[rp[[1]]==rp[[2]],
ReplacePart[c,{{rp[[1]]}}->2],
ReplacePart[c,(List/@rp)->1]],{n}]],
SparseArray[Table[ReplacePart[c,(List/@RandomSample[Range[p],2])->1],{n}]]]];

randomPedigree[{n_Integer,p_Integer},m_List,opts___Rule]:=Module[{k=Length[m]},
randomPedigree[{ConstantArray[n,k],ConstantArray[p,k]},m,opts]];

randomPedigree[{n:{___Integer},p:{___Integer}},m_List]:=Module[{demes=Range[Length[m]],ds,pars,c=ConstantArray[0,Total[p]],rp},
rp=partitionList[p];
ds=Flatten[MapThread[RandomChoice[#1->demes,{#2,2}]&,{m,n}],1];
pars=If[#[[1]]==#[[2]],RandomSample[rp[[#[[1]]]],2],{RandomChoice[rp[[#[[1]]]]],RandomChoice[rp[[#[[2]]]]]}]&/@ds;
SparseArray[ReplacePart[c,(List/@#)->1]&/@pars]];

Parents

This works, but relies on ArrayRules returning a sorted list of the parents, and on each individual having  1 or 2 parents. Not good coding...

Parents::usage="Parents[P] returns the parents of each individual in the pedigree.  P must be a SparseArray.";

Parents[P_SparseArray]:=Partition[#[[1,2]]&/@Drop[ArrayRules[P]/.{(x_->2)->Sequence[x->2,x->2 ]},-1],2];

TakeAway

TakeAway::usage="TakeAway[U,V]  gives the set U, minus elements V. Note that this is not the same as Complement[U,V] if U contains repeated elements.  ";

TakeAway[U_List,{}]:=U;
TakeAway[U_List,{v_}]:=Drop[U,Position[U,v][[1]]]/;Complement[{v},U]=={};
TakeAway[U_List,V_List]:=TakeAway[TakeAway[U,Take[V,1]],Drop[V,1]]/;Complement[V,U]=={};
TakeAway[U_List,W__List,V_List]:=TakeAway[TakeAway[U,W],V];

partitionList

partitionList::usage="partitionList[{(Subscript[a, 1]),…},{(Subscript[n, 1]),(Subscript[n, 2]),…}] partitions the list into sublists with {(Subscript[n, 1]),…} elements. partitionList[{(Subscript[n, 1]),(Subscript[n, 2]),…}] partitions {1,2,…}. partitionList[{(Subscript[a, 1]),…},{{(Subscript[n, 1, 1]),…},{(Subscript[n, 2, 1]),…},…} is also defined";

partitionList::badNumber="List `1` does not have the same length as the total of `2`";

partitionList[n:{___Integer}]:=partitionList[Range[Total[n]],n];
partitionList[a_List,n:{___Integer}]:=Module[{fl,nx},
If[Length[a]≠Total[n],Message[partitionList::badNumber,a,n]];
fl=FoldList[Plus,0,n];
nx=Transpose[{Drop[fl+1,-1],Drop[fl,1]}];
a[[Span@@#]]&/@nx];

partitionList[n:{{___Integer}...}]:=partitionList[Range[Total[n,2]],n];
partitionList[a_List,n:{{___Integer}...}]:=Module[{fpl},
fpl=partitionList[a,Flatten[n,1]];
partitionList[fpl,Length/@n]];

Discrete locus simulation

makeFamily

makeFamily::usage="makeFamily[(Subscript[ind, 1]),(Subscript[ind, 2]),k] makes a family of k offspring";

makeFamily[ind1:{{(0|1)...},{(0|1)...}},ind2:{{(0|1)...},{(0|1)...}},nkids_Integer]:=
Module[{mask,nloci=Length[ind1[[1]]]},
mask=RandomInteger[{0,1},{nkids,2,nloci}];
{(ind1[[1]](1-#[[1]])+ind1[[2]]#[[1]]),(ind2[[1]](1-#[[2]])+ind2[[2]]#[[2]])}&/@mask];

setαh

setαh::usage="setαh[(Subscript[n, loci]),γ] chooses α=1/(Sqrt[n]), and h=0 or 1 with probability (1/2)+(γ/(Sqrt[n])). This ensures that (Subscript[V, A]), (Subscript[V, D]), ℐ are all of order 1.";

setαh[n_Integer,γ_]:=Module[{k=Round[(n/2)+γ(Sqrt[n])]},
Transpose[{ConstantArray[(1/(Sqrt[n])),n],Join[ConstantArray[0,n-k],ConstantArray[1,k]]}]];

pop

pop::usage="pop[𝒫] lists the successive populations of diploid individuals for simulation 𝒫";

time

time::usage="time[𝒫] tracks the generation";

pedigree

pedigree::usage="pedigree[𝒫] stores a list of the full pedigrees associated with the symbol 𝒫.  This allows functions such as Fx to refer to the list without storing it.  It needs to be defined before use.  ** Here, a pedigree is represented as a sparse array of the full matrix - a different form from pedigreeList";

iterate

iterate::usage="iterate[𝒫,n,W] generates a pedigree and a set of n offspring, appending these to pedigreeList[𝒫], pop[𝒫].  iterate[𝒫,{(Subscript[n, 1]),…,(Subscript[n, T])},W] or iterate[𝒫,n,T,W] iterates over multiple generations. Selfing->True allows selfing. ";

Options[iterate]={Selfing->False};

iterate[𝒫_,n_Integer,W_,opts___Rule]:=Module[{ped, kids},
ped=selectPedigree[n,Last[pop[𝒫]],W,opts];
time[𝒫]+=1;
AppendTo[pedigree[𝒫],ped];
AppendTo[pedigreeList[𝒫],Parents[ped]];
AppendTo[pop[𝒫],makeKids[Last[pop[𝒫]],ped]]];

iterate[𝒫_,n_Integer,t_Integer,W_,opts___Rule]:=iterate[𝒫,ConstantArray[n,t],W,opts];
iterate[𝒫_,nl:{___Integer},W_,opts___Rule]:=Scan[iterate[𝒫,#,W,opts]&,nl];

initialise

initialise::usage="initialise[𝒫,N,(Subscript[n, loci])] sets up a simulation, 𝒫, with an initial population pop[𝒫] of N diploid individuals, each with (Subscript[n, loci]) loci, and allelic states labelled by ancestry (i.e., by 1…2N). Also sets up the pedigree list as {}. initialise[𝒫,{(Subscript[j, 1]),…}] initialises alleles {0,1} with (Subscript[j, i]) copies at each locus, 0≤j≤2ninds, in LE. initialise[𝒫,(Subscript[n, loci]),j] has the same initial # at all loci.";

initialise[𝒫_,pars__]:=Module[{},
time[𝒫]=0;pedigree[𝒫]={};pedigreeList[𝒫]={};pop[𝒫]={initPop[pars]}];

initPop

initPop::usage="initPop[N,(Subscript[n, loci])] sets up a population of N diploid individuals, each with (Subscript[n, loci]) loci, and allelic states labelled by ancestry (i.e., by 1…2N). initPop[N,{(Subscript[j, 1]),…}] initialises alleles {0,1} with (Subscript[j, i]) copies at each locus, 0≤j≤2ninds, in LE. initPop[N,(Subscript[n, loci]),j] has the same initial # at all loci.";

initPop[ninds_Integer,nloci_Integer]:=Module[{i},
Partition[Table[ConstantArray[i,nloci],{i,2ninds}],2]];

initPop[ninds_Integer,nloci_Integer,j_Integer]:=initPop[ninds,ConstantArray[j,nloci]];

initPop[ninds_Integer,jl:{___Integer}]:=Partition[Transpose[RandomSample[Join[ConstantArray[0,2ninds-#],ConstantArray[1,#]]]&/@jl],2];

selectPedigree

selectPedigree::usage="selectPedigree[n,parents,W] makes a random WF pedigree (but with no selfing), with n offspring; W[{{},{}}] is applied to diploid individuals to give the relative fitness. There must be at least two parents to avoid selfing. Returns a SparseArray. Selfing->True allows selfing, giving the WF model.";

Selfing::usage="Selfing is an option for randomPedigree and selectPedigree that allows selfing";

Options[selectPedigree]={Selfing->False};

selectPedigree[n_Integer,pars:{{{___Integer},{___Integer}}...},W_,opts___Rule]:=
Module[{wp=W/@pars,np=Length[pars],rp,c},
c=ConstantArray[0,np];
If[Selfing/.{opts}/.Options[selectPedigree],
SparseArray[Table[rp=RandomChoice[wp->Range[np],{2}];
If[rp[[1]]==rp[[2]],
ReplacePart[c,{{rp[[1]]}}->2],
ReplacePart[c,(List/@rp)->1]],{n}]],
SparseArray[Table[ReplacePart[c,(List/@RandomSample[wp->Range[np],2])->1],{n}]]]];

makeKids

makeKids::usage="makeKids[parents,P] generates offspring from the population, according to the pedigree P. Assumes no linkage.";

makeKids[pop:{{{___Integer},{___Integer}}...},ped:(_SparseArray|_List)]:=
Module[{pp,mask,pars,nkids=Length[ped],nloci=Length[pop[[1,1]]]},
pars=If[Max[#]==2,pp=Pick[pop,#,2];Join[pp,pp],Pick[pop,#,1]]&/@ped;
mask=RandomInteger[{0,1},{nkids,2,nloci}];
MapThread[(#1[[1]](1-#2)+#1[[2]]#2)&,{pars,mask},2]];

shuffle

shuffle::usage="shuffle[{{{0,1,…},{0,0,…}},…}] shuffles a diploid population, to simulate LE";

shuffle[pop_List]:=Partition[Transpose[RandomSample/@Transpose[Flatten[pop,1]]],2];

inflateMatrix

inflateMatrix::usage="inflateMatrix[M,k] generates matrix where each row is repeated k times";

inflateMatrix[m:(_SparseArray|_List),k_Integer]:=Flatten[Transpose[Partition[Flatten[ConstantArray[m,k],1],Length[m]]],1];

makePedigrees

makePedigrees::usage="makePedigrees[𝒫,(Subscript[n, inds]),T] sets up random neutral pedigrees, stored in pedigree[𝒫], and lists the pairs of parents in pedigreeList[𝒫]. Selfing->True allows selfing";

makePedigrees[𝒫_,n_,T_,opts___Rule]:=Module[{t},
pedigree[𝒫]=Table[randomPedigree[{ninds,ninds}],{t,T}];
pedigreeList[𝒫]=Parents/@pedigree[𝒫];];

selectKids

selectKids::usage="selectKids[pop,pedigree,k,f] generates k offspring genotypes, {(Subscript[X, 1]),…,(Subscript[X, k])}, for each pair of parents specified by the pedigree, and then selects the single offspring that has largest f[X]. selectKids[pop,pedigree,1] or selectKids[pop,pedigree] uses no selection.";

selectKids[pop_List,ped:(_SparseArray|_List),1]:=selectKids[pop,ped];
selectKids[pop_List,ped:(_SparseArray|_List),1,f_]:=selectKids[pop,ped];
selectKids[pop_List,ped:(_SparseArray|_List)]:=makeKids[pop,ped];
selectKids[pop_List,ped:(_SparseArray|_List),k_Integer,f_]:=Module[{kids},
kids=makeKids[pop,inflateMatrix[ped,k]];
Flatten[MaximalBy[#,f,1]&/@Partition[kids,k],1]];

evolve

evolve::usage="evolve[{𝒫,(SuperStar[𝒫])},k,f] iterates over multiple generations, using selectKids[] to apply within-family selection based on f[X], and the pedigree stored in pedigree[𝒫]. The initial population is taken from First[pop[𝒫]]. The results are appended to pop[(SuperStar[𝒫])]; thus, one could label replicates by successively calling . evolve[{𝒫,(SuperStar[𝒫])}] does not select.";

evolve[{𝒫_,𝒫r_}]:=Module[{ t,T=Length[pedigree[𝒫]]},
pop[𝒫r]={First[pop[𝒫]]};
Do[AppendTo[pop[𝒫r],selectKids[Last[pop[𝒫r]],pedigree[𝒫][[t]]]],{t,1,T}]];

evolve[{𝒫_,𝒫r_},k_Integer,f_]:=Module[{ t,T=Length[pedigree[𝒫]]},
pop[𝒫r]={First[pop[𝒫]]};
Do[AppendTo[pop[𝒫r],selectKids[Last[pop[𝒫r]],pedigree[𝒫][[t]],k,f]],{t,1,T}]];

randomJ

randomJ::usage="randomJ[n,k] chooses k integers on {1,n-1}, with probability proportional to 1/(j(n-j)). randomJ[n,p,F,k] samples k integers on {1,n-1}, from a large population with beta distributed allele frequencies.";

randomJ[n_Integer,k_Integer]:=Module[{rn=Range[1,n-1]},RandomChoice[(1/(rn(n-rn)))->rn,k]];

randomJ[n_Integer,p_,F_,k_Integer]:=Module[{rn=Range[1,n-1]},RandomChoice[BetaBinomial[rn,n,p,F]->rn,k]];

This samples integers from 1/(j(n-j)). An alternative is to sample from a Beta-binomial distribution:

( ({{n}, {j}}))((ℬ[j+α (OverBar[p]),n-j+α (OverBar[q])])/(ℬ[α (OverBar[p]),α (OverBar[q])]))   α=(1/F)-1

BetaBinomial

BetaBinomial::usage="BetaBinomial[j,n,(OverBar[p]),F] gives the PDF of a beta-binomial distribution";

BetaBinomial[j_,n_Integer,p_,F_]:=Module[{α=(1/F)-1},
Binomial[n,j]((Beta[j+α p,n-j+α(1-p)])/(Beta[α p,α(1-p)]))];

Components of mean and variance

z

z::usage="z[{{1,0,…},{0,1,…}},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] determines the trait value, with {-α,(2h-1)α,α} at each locus. z[{{{1,0,…},{0,1,…}},…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] applies to a list";

zOld[ind:{{(0|1)...},{(0|1)...}},αh:{{_,_}...}]:=
Total[MapThread[(#2[[1]](-1-#1+((#1)^2)+2#2[[2]]#1(2-#1)))&,{Total[ind],αh}]];(* slow *)

z[ind:{{(0|1)...},{(0|1)...}},αh:{{_,_}...}]:=Module[{α=αh[[All,1]],h=αh[[All,2]],x=Total[ind]},
(-1+x).α+(x(2-x)).(α(2h-1))];

z[inds:{{{(0|1)...},{(0|1)...}}...},αh:{{_,_}...}]:=Module[{α=αh[[All,1]],h=αh[[All,2]],xx=Total/@inds},
(-1+xx).α+(xx(2-xx)).(α(2h-1))];

z[{{0,1,1},{0,0,1}},{{(Subscript[α, 1]),(Subscript[h, 1])},{(Subscript[α, 2]),(Subscript[h, 2])},{(Subscript[α, 3]),(Subscript[h, 3])}}]//Simplify

{Timing[zOld[#,αh]&/@pop[sim][[-1]];],Timing[z[pop[sim][[-1]],αh];]}

z[pop[sim][[-1]],αh]==(zOld[#,αh]&/@pop[sim][[-1]])

zA

zA::usage="zA[{{1,0,…},{0,1,…}},{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] determines the additive component of the trait value, with {-α,(2h-1)α,α} at each locus; this is a(ζ+(SuperStar[ζ])), where a=2α(p(1-h)+q h). zA[{{{1,0,…},{0,1,…}},…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] applies to a list";

zA[ind:{{(0|1)...},{(0|1)...}},p_List,αh:{{_,_}...}]:=Module[{α=αh[[All,1]],h=αh[[All,2]],x=Total[ind]},
(x-2p).(2α(p(1-h)+(1-p)h))];

zA[inds:{{{(0|1)...},{(0|1)...}}...},p_List,αh:{{_,_}...}]:=Module[{α=αh[[All,1]],h=αh[[All,2]],xx=Total/@inds},
((#-2p)&/@xx).(2α(p(1-h)+(1-p)h))];

zD

zD::usage="zD[{{1,0,…},{0,1,…}},{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] determines the dominance component of the trait value, with {-α,(2h-1)α,α} at each locus; this is d((SuperStar[ζζ])), where a=2α(p(1-h)+q h). zD[{{{1,0,…},{0,1,…}},…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] applies to a list";

zD[ind:{{(0|1)...},{(0|1)...}},p_List,αh:{{_,_}...}]:=Module[{α=αh[[All,1]],h=αh[[All,2]],zz=(ind[[1]]-p)(ind[[2]]-p)},
zz.(-2α(2h-1))];

zD[inds:{{{(0|1)...},{(0|1)...}}...},p_List,αh:{{_,_}...}]:=Module[{α=αh[[All,1]],h=αh[[All,2]],zz=((#[[1]]-p)(#[[2]]-p))&/@inds},
zz.(-2α(2h-1))];

The difference between (Subscript[z, A])+(Subscript[z, D]) and the actual genotypic value is equal to the initial α(p-q)+2α (2h-1)p q - which is not the same as the initial mean.

z[pop[sim][[1]],αh]-(zA[pop[sim][[1]],pp,αh]+zD[pop[sim][[1]],pp,αh])

z[pop[sim][[-1]],αh]-(zA[pop[sim][[-1]],pp,αh]+zD[pop[sim][[-1]],pp,αh])

{α.(2pp-1)+(2α pp(1-pp)).(2h-1),Mean[z[pop[sim][[1]],αh]]}//N

zbar

zbar::usage="zbar[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives the expected trait mean";

zbar[p_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
α.(2p-1)+(2α p(1-p)).(2h-1)];

Va

Va::usage="Va[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives the additive variance. Va[{(Subscript[p, 1]),…},{(Subscript[p, 1]^ref),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] uses a different reference distribution.";

Va[p_List,αh_List]:=Va[p,p,αh];
Va[p_List,pr_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
8(p(1-p)).(((α(pr(1-h)+(1-pr)h)))^2)];

Vd

Vd::usage="Vd[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives the dominance variance. Vd[{(Subscript[p, 1]),…},{(Subscript[p, 1]^ref),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] uses a different reference distribution.";

Vd[p_List,αh_List]:=Vd[p,p,αh];
Vd[p_List,pr_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
(p(1-p)(p(1-p)+2(((p-pr))^2))).(((-2α(2h-1)))^2)];

Vad

Vad::usage="Vad[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives the covariance between A, D, which is zero if the allele frequencies equal the reference. Vad[{(Subscript[p, 1]),…},{(Subscript[p, 1]^ref),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] uses a different reference distribution.";

Vad[p_List,αh_List]:=0;
Vad[p_List,pr_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
-8(p(1-p)(p-pr)((α^2)(pr(1-h)+(1-pr)h))).((2h-1))];

Vg

Vg::usage="Vg[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives the total genotypic variance. Vg[{(Subscript[p, 1]),…},{(Subscript[p, 1]^ref),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] uses a different reference distribution.";

Vg[p_List,αh_List]:=Vg[p,p,αh];
Vg[p_List,pr_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
2(p(1-p)).((α^2)((((2h-1))^2)(1-2p(1-p))-2(2h-1)(2p-1)+1))];

Va[{p},{{α,h}}]+Vd[{p},{{α,h}}]-Vg[{p},{{α,h}}]//Simplify

Checking the algebra:

m={(q^2),2p q,(p^2)}.{-α,(2h-1)α,α};
v={(q^2),2p q,(p^2)}.((({-α,(2h-1)α,α}-m))^2);
8(α^2)p q(((p(1-h)+q  h))^2)+4(α^2)(p^2)(q^2)(((2h-1))^2)-v/.q->1-p//Simplify

ℐ

ℐ::usage="ℐ[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives the inbreeding depression. ℐ[{(Subscript[p, 1]),…},{(Subscript[p, 1]^ref),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] uses a different reference distribution.";

ℐ[p_List,αh_List]:=ℐ[p,p,αh];
ℐ[p_List,pr_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
(p(1-p)).(-2α(2h-1))];

ℐ2=(Subscript[V, D]) for biallelic loci

ℐ2::usage="ℐ2[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives ∑(Subscript[d, i]^2)((Subscript[p, i])(Subscript[q, i]))^2. ℐ2[{(Subscript[p, 1]),…},{(Subscript[p, 1]^ref),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] uses a different reference distribution.";

ℐ2[p_List,αh_List]:=ℐ2[p,p,αh];
ℐ2[p_List,pr_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
(((p(1-p)))^2).(((-2α(2h-1)))^2)];(* to fix *)

σADI

σADI::usage="σADI[p,αh] gives the covariance between A and D in inbred individuals";

σADI[p_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
2(p(1-p)(1-2p)).(-2α(2h-1)(2α(p(1-h)+(1-p)h)))];

σ2DI

σ2DI::usage="σ2DI[p,αh] gives the variance of inbreeding depression in inbred individuals";

σ2DI[p_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
(p(1-p)).((((-2α(2h-1)))^2))-4Vd[p,αh]];

zbarA

zbarA::usage="zbarA[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives the additive component of the mean. zbarA[{(Subscript[p, 1]),…},{(Subscript[p, 1]^ref),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] uses a different reference distribution.";

zbarA[p_List,αh_List]:=zbarA[p,p,αh];
zbarA[p_List,pr_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
2α.(p(1-h)+(1-p)h)];

zbarD

zbarD::usage="zbarD[{(Subscript[p, 1]),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] gives the dominance component of the mean. zbarD[{(Subscript[p, 1]),…},{(Subscript[p, 1]^ref),…},{{(Subscript[α, 1]),(Subscript[h, 1])},…}] uses a different reference distribution.";

zbarD[p_List,αh_List]:=zbarD[p,p,αh];
zbarD[p_List,pr_List,αh_List]:=Module[{α=αh[[All,1]],h=αh[[All,2]]},
-2α.(2h-1)];

({{X, 0, 1, 2}, {ζ+(SuperStar[ζ]), -2p, q-p, 2q}, {ζ+(SuperStar[ζ]), (p^2), -p q, (q^2)}, {z, -α, (2h-1)α, α}})         ({{𝔼[z]=α(p-q)+2p q(2h-1)}, {𝔼[z(ζ+(SuperStar[ζ]))]=2a p q=4α p q(p(1-h)+q h)}, {𝔼[z((SuperStar[ζζ]))]=d (p^2)(q^2)=-2α(2h-1)(p^2)(q^2)}})   hence     ({{a=2α (p(1-h)+q h), (Subscript[V, A])=8(∑_i)(Subscript[α, i]^2)(Subscript[p, i])(((Subscript[q, i])((Subscript[p, i])(1-(Subscript[h, i]))+(Subscript[q , i])(Subscript[h, i])))^2)}, {d=-2α(2h-1), (Subscript[V, D])=4(∑_i)(Subscript[α, i]^2)(Subscript[p, i]^2)(((Subscript[q, i]^2)(2(Subscript[h, i])-1))^2)}, { , (Subscript[V, A])+(Subscript[V, D])=2(∑_i)(Subscript[α, i]^2)(Subscript[p, i])(Subscript[q, i])((((2(Subscript[h, i])-1))^2)(1-2(Subscript[p, i])(Subscript[q, i]))-2(2(Subscript[h, i])-1)((Subscript[p, i])-(Subscript[q, i]))+1)}, { , ℐ=(∑_i)(Subscript[d, i])(Subscript[p, i])(Subscript[q, i])=-2(∑_i)(Subscript[α, i])(Subscript[p, i])(Subscript[q, i])(2(Subscript[h, i])-1)}})

makeDiploidEffects

Not clear that this is a sensible approach

makeDiploidEffects::usage="makeDiploidEffects[n,j,h] generates an n×n matrix of diploid genotypic values, assuming n-j - alleles, j + alleles, such that (Subscript[g, 0, 0])=-1, (Subscript[g, 0, 1])=2h -1, (Subscript[g, 1, 1])=1. makeDiploidEffects[n,{(Subscript[h, 1]),…},{(Subscript[α, 1]),…}] has effects (Subscript[g, i, i])=(Subscript[α, i]) in homozygotes and dominance determined by the ordering of the (Subscript[h, i]). makeDiploidEffects[n,(Subscript[v, a]),(Subscript[v, d]),θ] draws additive effects from a Gaussian, variance (Subscript[v, a]), and a random dominance effect from a Gaussian, variance (Subscript[v, d]), plus a mean (Subscript[d, i, i])=θ (Subscript[α, i]) ";

makeDiploidEffects[n_Integer,j_Integer,h_]:=ArrayFlatten[{{ConstantArray[1,{j,j}],ConstantArray[2h-1,{j,n-j}]},{ConstantArray[2h-1,{n-j,j}],ConstantArray[-1,{n-j,n-j}]}}]

makeDiploidEffects[n_Integer,h_List,α_List]:=Module[{i,j},
Table[If[h[[i]]>h[[j]],α[[i]],α[[j]]],{i,n},{j,n}]];

makeDiploidEffects[n_Integer,va_,vd_,θ_]:=Module[{a},
a=RandomReal[NormalDistribution[0,(Sqrt[va])],n];
Outer[Plus,a,a]+θ DiagonalMatrix[a]+RandomReal[NormalDistribution[0,(Sqrt[vd])],{n,n}]];

makeDiploidEffects[10,1,1,0.5]//MatrixForm

Replicate simulations with selection

makeSel

makeSel::usage="makeSel[{𝒩,𝒫},M,(Subscript[k, s]),(Subscript[n, reps]),(Subscript[n, inds]),T] sets up (Subscript[n, reps]) replicate simulations, based on pedigree 𝒩 and initial frequencies prM[M] and effects αhM[M]; truncation selection is 1/(Subscript[k, s]).  These generate data labelled  by 𝒫[M,(Subscript[k, s]),(Subscript[n, reps])][rep], and store zAlist[], zDList[], zbAlist[], zbDlist[], vAlist[], vDlist[], vADlist[],vAGlist[], vDGlist[], vADGlist[]. Does not store pop[] or pList[]. ";

makeSel[{𝒩_,𝒫_},M_Integer,ks_Integer,nreps_Integer,ninds_Integer,T_Integer]:=Module[{jrr,𝒳},
Do[jrr=RandomInteger[BinomialDistribution[2ninds,#]]&/@prM[M];
𝒳=𝒫[M,ks,nreps][r];
pop[𝒩]={initPop[ninds,jrr]};
evolve[{𝒩,𝒳},ks,z[#,αhM[M]]&]; 
 pList[𝒳] = N[Mean[Flatten[#, 1]]] & /@ pop[𝒳];
zAlist[𝒳] = zA[#, prM[M], αhM[M]] & /@ pop[𝒳];
zDlist[𝒳] = zD[#, prM[M], αhM[M]] & /@ pop[𝒳];
zbAlist[𝒳]=Mean/@zAlist[𝒳];zbDlist[𝒳]=Mean/@zDlist[𝒳];
vAlist[𝒳]=Variance/@zAlist[𝒳];vAGlist[𝒳]=Va[#,prM[M], αhM[M]]&/@pList[𝒳];
vDlist[𝒳]=Variance/@zDlist[𝒳];vDGlist[𝒳]=Vd[#,prM[M], αhM[M]]&/@pList[𝒳];
vADlist[𝒳]=MapThread[Covariance,{zAlist[𝒳],zDlist[𝒳]}];vADGlist[𝒳]=Vad[#,prM[M], αhM[M]]&/@pList[𝒳];pop[𝒳] =.;pList[𝒳] =.;,
 {r, nreps}];]

zAlist, zDlist,zbGlist, zbAlist, zbDlist, vGlist,vAlist, vDlist, vADlist,vGGlist, vAGlist, vDGlist, vADGlist

zAlist::usage="zAlist[𝒫] stores A for each individual";
zDlist::usage="zDlist[𝒫] stores D for each individual";
zbAlist::usage="zbAlist[𝒫] stores mean A for each generation";
zbDlist::usage="zbDlist[𝒫] stores mean D for each generation";

vAlist::usage="vAlist[𝒫] stores var(A) for each generation";
vDlist::usage="vDlist[𝒫] stores var(D) for each generation";
vADlist::usage="vADlist[𝒫] stores cov(A,D) for each generation";
vAGlist::usage="vAGlist[𝒫] stores the genic component of var(A) for each generation";
vDGlist::usage="vDGlist[𝒫] stores the genic component of var(D) for each generation";
vADGlist::usage="vADGlist[𝒫] stores the genic component of cov(A,D) for each generation";

zbGlist[𝒫_]:=zbAlist[𝒫]+zbDlist[𝒫];
vGlist[𝒫_]:=vAlist[𝒫]+vDlist[𝒫]+2vADlist[𝒫];
vGGlist[𝒫_]:=vAGlist[𝒫]+vDGlist[𝒫]+2vADGlist[𝒫];

fPlot

fPlot::usage="fPlot[𝒫,{(Subscript[f, 1]),…},{(Subscript[c, 1]),…},(Subscript[n, reps]),T,{(Subscript[y, 0]),(Subscript[y, 1])},(Subscript[y, label])] plots the specified functions (Subscript[f, i])[𝒫[r]] in the correspondng styles, averaging over replicates. fPlot[𝒫,{(Subscript[f, 1]),…},{(Subscript[c, 1]),…},{1,2,…},T,{(Subscript[y, 0]),(Subscript[y, 1])},(Subscript[y, label])] plots the specified replicates. Can add options.";

fPlot[𝒫_,f_List,c_List,nreps_Integer,T_Integer,{y0_,y1_},ylab_,opts___Rule]:=Module[{r,fl},
fl=(#[𝒫[r]]&)/@f;
Show[MapThread[ListLinePlot[#1,PlotStyle->#2,PlotRange->{{0,T+1},{y0,y1}},BaseStyle->{FontFamily->"Times",FontSize->"14"},AxesLabel->{"time",ylab},opts]&,
{Mean[Table[fl,{r,nreps}]],c}]]];

fPlot[𝒫_,f_List,c_List,reps_List,T_Integer,{y0_,y1_},ylab_,opts___Rule]:=Module[{r,fl,j},
Show[Table[fl=(#[𝒫[reps[[r]]]]&)/@f;
MapThread[ListLinePlot[#1,PlotStyle->#2,PlotRange->{{0,T+1},{y0,y1}},BaseStyle->{FontFamily->"Times",FontSize->"14"},AxesLabel->{"time",ylab},opts]&,{fl,c}],{r,Length[reps]}]]];