Nothing
## Ranking method
ranking.strains=function(DGobject, bw, nb.mcsimul, plots=FALSE,
kernel.type="Quadratic"){
if(!(kernel.type%in%c("Linear","Quadratic","Power3","Power4"))){
kernel.type="Quadratic"
print("WARNING: specified kernel.type is unknown. Default Quadratic kernel is used.")
}
geneticData=DGobject@genetic
x=DGobject@demographic[,1:2]
Z=DGobject@demographic[,3]
nbStrains=ncol(geneticData)-2
N=nrow(x)
## Estimation of strain proportions
pSEstimes=.estimation.pS(geneticData,x,bw,kernel.type)
## Estimation of regression coefficients
regression <- lm(Z ~ -1 + pSEstimes)
z=regression$coef
## Assessment of the significance of differences in the coefficients
zStar=matrix(0,nb.mcsimul,length(z))
for(i in 1:nb.mcsimul){
permut=(1:N)[!is.na(pSEstimes[,1])]
pSEstimesStar=pSEstimes
pSEstimesStar[permut,]=pSEstimesStar[sample(permut),]
regressionStar=lm(Z ~ -1 + pSEstimesStar)
zStar[i,]=regressionStar$coefficients
}
## Computation of pairwise t-tests
pvals=NULL
pvals.name=NULL
for(i in 1:(ncol(zStar)-1)){
for(j in (i+1):ncol(zStar)){
pvals=c(pvals,mean( zStar[,j]-zStar[,i] >= z[j]-z[i] ))
pvals.name=c(pvals.name,paste("(",j,"-",i,")>0",sep=""))
pvals=c(pvals,mean( zStar[,j]-zStar[,i] <= z[j]-z[i] ))
pvals.name=c(pvals.name,paste("(",j,"-",i,")<0",sep=""))
pvals=c(pvals,mean( abs(zStar[,j]-zStar[,i]) >= abs(z[j]-z[i]) ))
pvals.name=c(pvals.name,paste("(",j,"-",i,")neq0",sep=""))
}
}
names(pvals)=pvals.name
if(plots){
## Production of plots (ranking plots)
nrow.plot=ceiling((2+nbStrains)/3)
par(mfrow=c(nrow.plot,3),mar=c(2,2,1,0.2))
plot(x,cex=(Z-min(Z))/(max(Z)-min(Z))*3,pch=19,asp=1,axes=F,xlab="",ylab="",
main="Growth variable")
plot(x,cex=1,pch=1,asp=1,axes=F,xlab="",ylab="",main="Sampling sites")
points(geneticData[,1:2],pch=19)
for(j in 1:nbStrains){
zseq=seq(min(c(z,zStar),na.rm=TRUE),max(c(z,zStar),na.rm=TRUE),l=20)
hist(zStar[,j],xlim=range(zseq),breaks=zseq,xlab="",ylab="",main=paste("Coef. strain ",j))
abline(v=z[j],lwd=2,lty="dashed",col=2)
}
}
names(z)=paste("strain",1:nbStrains,sep="")
colnames(zStar)=paste("strain",1:nbStrains,sep="")
return(list(permutation.estimates=zStar,estimates=z,p.values=pvals))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.