#######################################################
### Functions
.rateClonemates1G=function(counts1){
return(sum(counts1/sum(counts1)*pmax(0,counts1-1)/(sum(counts1)-1)))
}
.rateClonemates2G=function(counts1,counts2){
return(sum(counts1/sum(counts1)*counts2/sum(counts2)))
}
.count2profil=function(R){
out=NULL
for(i in 1:length(R)){ out=c(out,rep(i,R[i])) }
return(out)
}
.profil2count=function(Z){
out=NULL
for(i in 1:length(unique(Z))){ out=c(out,sum(Z==unique(Z)[i])) }
return(out)
}
.cycle=function(R0,s){
Z0=.count2profil(R0)
N=length(Z0)
S=rbinom(1,N,s)
C=N-S
Z1C=sample(Z0,C,replace=T)
R1=NULL
for(i in 1:length(R0)){ R1=c(R1,sum(Z1C==i)) }
return(cbind(c(R0,rep(0,S)),c(R1,rep(1,S))))
}
.cycleMult=function(R0,s,n1,n2){
Ri=R0
if(n1>0){
for(i in 1:n1){
Ri=.cycle(Ri,0)[,2]
}
}
Rn1plus1=.cycle(Ri,s)[,2]
R0=c(R0,rep(0,length(Rn1plus1)-length(R0)))
Ri=Rn1plus1
if(n2>0){
for(i in 1:n2){
Ri=.cycle(Ri,0)[,2]
}
}
return(cbind(R0,Ri))
}
.logit=function(p){ return(log(p/(1-p))) }
.logitinv=function(x){ return(exp(x)/(1+exp(x))) }
## Contrast
.critMult=function(f0raw,f1raw,p0,p1,g01,n1,n2,logN,logits){
N=round(exp(logN))
s=.logitinv(logits)
f0=N*f0raw/(N-1)-p0/(N-1)
Fobs=N*f1raw/(N-1)-p1/(N-1)
Gobs=g01
a=1/N
b=1-1/N
G=(1-s)*(a+b*f0)
fprime=(1-s)^2*(a*sum(b^(0:n1))+b^(n1+1)*f0)
if(n2>0){
fprime=a*sum(b^(0:(n2-1)))+b^n2*fprime
}
return(abs(fprime-Fobs)+abs(G-Gobs))
}
## Optimization function
.optim.Ne.s=function(f0raw,f1raw,p0,p1,g01,n1,n2,initial,maxit){
h=function(par){
.critMult(f0raw,f1raw,p0,p1,g01,n1=n1,n2=n2,par[1],par[2])
}
oo=NULL
oo$value=Inf
for(i in 1:nrow(initial)){
ooNew=optim(c(log(initial[i,1]),.logit(initial[i,2])),h,
control=list(maxit=maxit))
if(ooNew$value<oo$value){ oo=ooNew }
}
oo$partransf=c(round(exp(oo$par[1])),.logitinv(oo$par[2]))
return(oo)
}
## Simulation
simul=function(start.counts,param,nb1ClonalCycles,nb2ClonalCycles,nb1Sample,nb2Sample,
nbPrelCycles){
Ne=param[1]
s=param[2]
## initial conditions
if(sum(start.counts)<Ne){
start.counts=start.counts[start.counts>0]
start.countsstar=c(sum(start.counts==1),start.counts[start.counts>1])
R0=as.numeric(rmultinom(1,Ne,start.countsstar))
R0=c(R0[-1],rep(1,R0[1]))
}
if(sum(start.counts)==Ne){
R0=start.counts
}
if(sum(start.counts)>Ne){
R0=as.numeric(rmultinom(1,Ne,start.counts))
}
## R1 simulation
temp=.cycleMult(R0,s,nb1ClonalCycles,nb2ClonalCycles)
if(nbPrelCycles>0){
for(i in 1:nbPrelCycles){
temp=.cycleMult(temp[temp[,2]>0,2],s,nb1ClonalCycles,nb2ClonalCycles)
}
}
R0=temp[,1]
R1=temp[,2]
## sampling
z0=sample(.count2profil(R0),nb1Sample,replace=TRUE)
z1=sample(.count2profil(R1),nb2Sample,replace=TRUE)
r0=NULL
r1=NULL
for(i in sort(unique(c(z0,z1)))){
r0=c(r0,sum(z0==i))
r1=c(r1,sum(z1==i))
}
## statistics
f0raw=.rateClonemates1G(r0)
f1raw=.rateClonemates1G(r1)
p0=sum(r0/nb1Sample)
p1=sum(r1/nb2Sample)
g01=.rateClonemates2G(r0,r1)
#f0raw=(.rateClonemates1G(R0)+1/(Ne-1))*(Ne-1)/Ne
#f1raw=(.rateClonemates1G(R1)+1/(Ne-1))*(Ne-1)/Ne
#p0=sum(r0/nb1Sample)
#p1=sum(r1/nb2Sample)
#g01=.rateClonemates2G(R0,R1)
population=temp[apply(temp,1,function(x) sum(x)>0),]
colnames(population)=c("t0","t1")
samples=cbind(r0,r1)
colnames(samples)=c("t0","t1")
return(list(population=population,samples=samples,
statistics=c(f0raw,f1raw,p0,p1,g01)))
}
## parametric bootstrap
.parboot.parfix=function(r00,nech1,nech2,param,n1,n2,nbB,nbPrelCycles,param0,maxit,trace){
temp0=NULL
for(k in 1:nbB){
if(k/100==round(k/100) & trace){ print(paste("Iteration in parametric bootstrap:",k),quote=F) }
simk=simul(r00,param,n1,n2,nech1,nech2,nbPrelCycles)
f0raw=simk$statistics[1]
f1raw=simk$statistics[2]
p0=simk$statistics[3]
p1=simk$statistics[4]
g01=simk$statistics[5]
## parameter estimation
oo=.optim.Ne.s(f0raw,f1raw,p0,p1,g01,n1,n2,param0,maxit)
temp0=rbind(temp0,c(param,oo$partransf))
}
return(temp0)
}
.estim.CI=function(param,paramB,confidenceLevel){
## without bias correction
endpoints=c((1-confidenceLevel)/2,(1+confidenceLevel)/2)
est1=cbind(param,t(apply(paramB[,3:4],2,quantile,endpoints)))
## with bias correction
bias=c(mean(param[1]-paramB[,3]),mean(param[2]-paramB[,4]))
sdhat=apply(paramB[,3:4],2,sd)
param2=param+bias
est2=cbind(param2,rbind(qnorm(endpoints,param[1]+bias[1],sdhat[1]),
qnorm(endpoints,param[2]+bias[2],sdhat[2])))
est=cbind(est1,est2)
est[1,]=round(est[1,])
est[2,]=round(est[2,],digits=4)
rownames(est)=c("Ne","s")
colnames(est)=c("Estim","IC1","IC2","EstimBiasCorr","IC1BiasCorr",
"IC2BiasCorr")
return(est)
}
.estim.CItransf=function(param,paramB,confidenceLevel){
## without bias correction
endpoints=c((1-confidenceLevel)/2,(1+confidenceLevel)/2)
est1=cbind(param,t(apply(paramB[,3:4],2,quantile,endpoints)))
## with bias correction
##bias=c(mean(param[1]-paramB[,3]),mean(param[2]-paramB[,4]))
##sdhat=apply(paramB[,3:4],2,sd)
##param2=param+bias
##est2=cbind(param2,rbind(qnorm(endpoints,param[1]+bias[1],sdhat[1]),
## qnorm(endpoints,param[2]+bias[2],sdhat[2])))
## with bias correction
## after log transf for Ne and logit transformation for s
lparam=c(log(param[1]),log(param[2]/(1-param[2])))
lparamB=cbind(log(paramB[,3]),log(paramB[,4]/(1-paramB[,4])))
bias=c(mean(lparam[1]-lparamB[,1]),mean(lparam[2]-lparamB[,2]))
sdhat=apply(lparamB,2,sd)
lparam3=lparam+bias
browser()
est3=cbind(lparam3,rbind(qnorm(endpoints,lparam[1]+bias[1],sdhat[1]),
qnorm(endpoints,lparam[2]+bias[2],sdhat[2])))
est3=rbind(exp(est3[1,]),exp(est3[2,])/(1+exp(est3[2,])))
est=cbind(est1,est3)
est[1,]=round(est[1,])
est[2,]=round(est[2,],digits=4)
rownames(est)=c("Ne","s")
colnames(est)=c("Estim","IC1","IC2","EstimBiasCorr",
"IC1BiasCorr","IC2BiasCorr")
est
}
.coverage=function(param,nbCoverage,confidenceLevel,nbBootstrap,
r00,nech0,nech1,n1,n2,param0,maxit,trace){
cover=cbind(matrix(0,2,2),NA)
confint=NULL
colnames(cover)=c("Cover","CoverBiasCorr","nbIteration")
rownames(cover)=c("Ne","s")
for(k in 1:nbCoverage){
if(k/10==round(k/10) & trace){
print(paste("Iteration in coverage assessment:",k),quote=F)
cover[1,3]=k^2
}
oostar=.parboot.parfix(r00,nech0,nech1,param,
n1=n1,n2=n2,1,0,param0,maxit,trace=FALSE)
PBstar=.parboot.parfix(r00,nech0,nech1,oostar[3:4],
n1=n1,n2=n2,nbBootstrap,0,param0,maxit,trace=FALSE)
est=.estim.CI(oostar[3:4],PBstar,confidenceLevel)
cover[1,1]=cover[1,1]+((param[1]>=est[1,2]) & (param[1]<=est[1,3]))
cover[2,1]=cover[2,1]+((param[2]>=est[2,2]) & (param[2]<=est[2,3]))
cover[1,2]=cover[1,2]+((param[1]>=est[1,5]) & (param[1]<=est[1,6]))
cover[2,2]=cover[2,2]+((param[2]>=est[2,5]) & (param[2]<=est[2,6]))
confint=rbind(confint,c(est[1,2:3],est[2,2:3],est[1,5:6],est[2,5:6]))
}
cover[,1:2]=cover[,1:2]/nbCoverage
cover[1,3]=nbCoverage
return(cover)
}
###################################################
ncase=function(counts,nb1ClonalCycles,nb2ClonalCycles,start.param,maxit=500,
nbBootstrap=0,conf.level=0.95,nbCoverage=0,trace=TRUE){
counts1=counts[,1]
counts2=counts[,2]
## Computation of clonemate rates
f=.rateClonemates1G(counts1)
fprime=.rateClonemates1G(counts2)
p=sum(counts1/sum(counts1))
pprime=sum(counts2/sum(counts2))
g=.rateClonemates2G(counts1,counts2)
## Estimation
estim1=.optim.Ne.s(f,fprime,p,pprime,g,nb1ClonalCycles,nb2ClonalCycles,start.param,
maxit)
param1=estim1$partransf
est=param1
names(est)=c("Ne","s")
if(nbBootstrap>0){
print("PARAMETRIC BOOTSTRAP: BIAS CORRECTION AND CONFIDENCE INTERVALS (may be time consuming)",
quote=F)
estimB=.parboot.parfix(counts1,sum(counts1),sum(counts2),param1,
nb1ClonalCycles,nb2ClonalCycles,nbBootstrap,nbPrelCycles=0,start.param,maxit,
trace)
est=.estim.CI(param1,estimB,conf.level)
param2=est[,4]
if(nbCoverage>0){
print("ASSESSMENT OF THE COVERAGE OF CONFIDENCE INTERVALS (may be time consuming)",quote=F)
cover1=.coverage(param1,nbCoverage,conf.level,nbBootstrap,counts1,
sum(counts1),sum(counts2),nb1ClonalCycles,nb2ClonalCycles,start.param,
maxit,trace)
return(list(bootstrap.estim=estimB[,3:4],cover=cover1,
optim.output=estim1,estimates=est))
} else {
return(list(bootstrap.estim=estimB[,3:4],optim.output=estim1,estimates=est))
}
} else {
return(est)
}
}
ncase4simul=function(param,nb1ClonalCycles,nb2ClonalCycles,nb1Sample,nb2Sample,
start.counts,nbPrelCycles,nbSimul,start.param,maxit=500,trace=TRUE){
print("SIMULATIONS UNDER FIXED PARAMETERS AND RAW ESTIMATIONS (may be time consuming)",quote=F)
estimB=.parboot.parfix(start.counts,nb1Sample,nb2Sample,param,
nb1ClonalCycles,nb2ClonalCycles,nbB=nbSimul,nbPrelCycles,start.param,maxit,trace)
colnames(estimB)=c("Ne","s","Ne.hat","s.hat")
return(estimB)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.