Nothing
sim.several.s <- function
### Do several simulations each with only 1 s value.
(s,
### Vector of s values.
...
### Other arguments to sim.drift.selection.
){
mlply(data.frame(s),function(s)sim.drift.selection(s=s,...))
### List of selection simulations.
}
nicholsonppp.list <- function
### Fit a model for each simulation.
(sims
### List of simulations from nicholsonppp.
){
lapply(sims,function(L)nicholsonppp(L$sim[,,L$p$gen]))
### List of nicholson model fits.
}
ppp.df <- function
### Convert simulation and model fit lists to a data frame that can be
### used for plotting diagnostics for the ppp-value classifier.
(sims,
### List of simulations.
models
### List of nicholson model fits.
){
res.df <- ldply(1:length(models),
function(i)data.frame(ppp=models[[i]]$ppp,
type=sims[[i]]$s$type,
s=max(sims[[i]]$s$s)))
res.df$s <- factor(res.df$s,labels=format(unique(res.df$s),digits=2))
attr(res.df,"N") <- sims[[1]]$p$n.loc
res.df
### Data frame suitable for plotting densities or sensitivity.
}
dens.several.s <- function
### Plot density curves for ppp values for each selection type, to
### visualize the extent to which PPP-values can be used to indicate
### selection state in a simulation.
(df,
### Data frame with columns ppp type s to be plotted.
f=~ppp|s,
### Lattice plot formula for densityplot.
main="Small PPP-values indicate strong positive selection",
...
### To be passed to densityplot.
){
layout <- if(missing(f))c(1,nlevels(df$s)) else NULL
p <- densityplot(f,df,groups=type,
layout=layout,n=500,
strip=strip.custom(
strip.levels=c(TRUE,TRUE),strip.names=TRUE),
main=main,
...)
direct.label(p)
### The lattice plot.
}
hilite.best <- function
### panel.groups function for highlighting the best points by drawing
### grey lines and labeling the actual values.
(x,
y
){
panel.abline(v=x,col="grey")
panel.abline(h=y,col="grey")
first.x <- x[1]
grid.text(round(first.x,2),unit(first.x,"native"),1,
rot=90,just=c("right","bottom"),gp=gpar(col="grey"))
if(length(x)>1){
last.x <- x[length(x)]
grid.text(round(last.x,2),unit(last.x,"native"),1,
rot=90,just=c("right","top"),gp=gpar(col="grey"))
}
grid.text(round(y,2),1,unit(y,"native"),
just=c("right","bottom"),gp=gpar(col="grey"))
}
classify.loci <- function
### Classify loci into selection or not.
(res.df,
### Result of ppp.df.
xmax=0.6,
ylim=c(-2,ymax),
cutoff=seq(0,xmax,l=200),
### Vector of cutoff values to use for the classifier.
splitby=c("s")
){
classify.group <- function(sdf,cutoff=NULL){
i <- which(sdf$ppp<xmax)
if(is.null(cutoff))cutoff <- c(0,sort(sdf$ppp[c(i,max(i)+1)]))
classify1 <- function(cutoff){
summarise(transform(sdf,guess=ppp<cutoff,positive=type=="positive"),
true.positive=sum(guess==TRUE&positive==TRUE),
false.positive=sum(guess==TRUE&positive==FALSE),
true.negative=sum(guess==FALSE&positive==FALSE),
false.negative=sum(guess==FALSE&positive==TRUE),
sensitivity=sum(guess==TRUE&positive==TRUE)/sum(positive),
specificity=1-sum(guess==TRUE&positive==FALSE)/sum(!positive))
}
mdply(data.frame(cutoff),classify1)
}
if("desc"%in%names(res.df))splitby <- c(splitby,"desc")
ddply(res.df,splitby,classify.group,cutoff)
}
cutoff.plot <- function
### Plot prediction counts against cutoff values.
(cl,
### Data frame from classify.loci.
ylim=c(-5,21),
### Limits for y axis.
xlim=c(-0.15,0.6),
### Limits for x axis.
dens=NULL,
### Optional density data to plot.
f=percent~cutoff|s,
### Plot formula for xyplot.
...
### args for xyplot.
){
N <- sum(cl[1,3:6])
cl2 <- transform(cl,
correct=true.positive+true.negative,
incorrect=false.positive+false.negative)
molt <- transform(melt(cl2,id=1:2),percent=value/N*100)
minrisk.panel <- function(x,y,subscripts,group.number,...){
panel.xyplot(x=x,y=y,group.number=group.number,...)
if(group.number==4){
best.y <- min(y)
if(y[1]!=best.y){
best.x <- x[which(y==min(y))]
hilite.best(best.x,best.y)
}
}
}
## cool but useless:
##dl(xyplot,molt,value~cutoff,variable,type='l')
ss <- subset(molt,variable%in%c("true.positive","incorrect","false.positive","false.negative"))
ss$variable <- factor(ss$variable)
p2 <- xyplot(f,ss,
groups=variable,
panel=function(subscripts,groups,...){
panel.superpose(subscripts=subscripts,groups=groups,...)
if(!is.null(dens)){
d <- subset(dens,s==ss[subscripts,"s"][1])
panel.superpose.dl(d$ppp,NULL,1:nrow(d),d$type,
panel.groups="panel.rug",
end=0.05,
col=selection.colors.default)
}
},
panel.groups=minrisk.panel,
type='l',layout=c(1,nlevels(molt$s)),
xlim=xlim,
ylim=ylim,
ylab="Percent of loci in simulation",
xlab="Cutoff for PPP-value decision rule",
main="Prediction rates change with PPP-value cutoffs",
par.settings=list(superpose.line=list(col=c("brown","orange","violet","darkblue"))),
...)
direct.label(p2,method=function(...){
d <- first.points(...)
d$vjust <- c(1,0)
d
})
### The lattice plot.
}
roc.loci <- function
### Plot ROC curves of loci classification for several s values.
(cl,
### Result of classify.loci.
f=sensitivity~1-specificity|desc,
### Plot formula for xyplot.
...
### Arguments for xyplot.
){
roc <- xyplot(f,cl,
type='l',
panel=function(...){
panel.abline(0,1,col="grey")
panel.xyplot(...)
},
main="ROCs vary with selection strength",
...)
direct.label(roc,method=function(...){
x <- data.frame(perpendicular.lines(...),hjust=0)
x[2,c("x","y")] <- c(0.6,0.8)
x
})
### The lattice plot.
}
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.