varCKMR is a pseudo R package to showcase a method of estimating the variance in
the distribution of number of offspring. It does this using within birth-cohort
sibling comparisons.
It contains an individual based simulation that installs as the R package varCKMR along with
some helper functions.
The TMB models used in the analysis are in the models folder with the ssmodelAdj.cpp
being the model using within-cohort sibling comparisons and the ssmodelNoSC.cpp model
omitting them.
The scripts folder contains the scripts to recreate the simulations used in the analysis. populations
is where the populations used in the analysis would be kept althought they are ommitted here for space reasons.
sims contains all the samples from each population but are again ommitted for space reasons.
The org file varCKMR.org contains all the files to make all the files in this folder.
The code here uses a more complicated derivation to break up \(σ^2tot\) (and \(Var(R)\)) then what is given in the paper:
\begin{multline}
\label{muteqlaw2}
σ^2tot = Var(X) = ∑a=0^A Var(X|D = a)Δ_a + ∑a=0^A E[X| D = a]^2(1-Δ_a)Δ_a-
2∑a=1^A∑b=0a-1 E[X|D = a] Δ_a E[X|D = b] Δ_b.
\end{multline}
This form assumes that X is also mutually exclusive (which age of death and age of parent are) but if that is the case then the form above and the form given in the paper will agree.
##See README
##here is a small example using the sim parameters from one of the populations
surv = c(0.215,0.280)
death_probs = c(1-surv[1],surv[1]*(1-surv[2]),surv[1]*surv[2])
fecundity = 2*c(0.380,1.518,3.416)
theta = 0.1
EXD = cumsum(fecundity)
V1 = fecundity+fecundity^2/theta
varXD = cumsum(V1)
oldform <- function(varXD,EXD,Delta){
firstsum <- sum(varXD*Delta)
secondsum <- sum(EXD^2*(1-Delta)*(Delta))
thirdsum <- 0
for(i in 2:length(EXD)){
for(j in 1:(i-1)){
thirdsum = thirdsum + EXD[i]*Delta[i]*EXD[j]*Delta[j]
}
}
total <- firstsum + secondsum -2*thirdsum
total
}
newform <- function(varXD,EXD,Delta){
mu_tot = sum(EXD*Delta)
firstsum = sum(varXD*Delta)
secondsum = sum(EXD^2*Delta)
total = firstsum+secondsum-mu_tot^2
total
}
oldv = oldform(varXD,EXD,death_probs)
newv = newform(varXD,EXD,death_probs)
print(paste0("old: ",round(oldv,3)," new: ", round(newv,3)))old: 61.673 new: 61.673