etc/sensitivity.R

In nicholsonppp: Evolution simulation and classification

source("../R/sim.R")
source("../R/plot.R")
source("../R/sensitivity.R")
source("../R/model.R")
library(ggplot2)
library(lattice)
library(ff)
library(latticedl)
dyn.load("../src/0mers_twist.so")
## initial study of which s values show a good range of selection
## behavior
s <- 10^c(-3,-2,-1.5,-1,0)
sim <- sim.drift.selection(loci=100,generations=100,s=s,p.neutral=0.1)
df <- sim2df(sim)
sim.summary.plot(df)



sims <- sim.several.s(s,loci=100,generations=100)
models <- nicholsonppp.list(sims)
##save(sims,models,file="sims.models.Rdata")
load(file="sims.models.Rdata")



sims.neu <- sim.several.s(s,loci=100,generations=100,
                          p.neutral=0.99,array.fun=ff)
models.neu <- nicholsonppp.list(sims.neu)
##save(sims.neu,models.neu,file="sims.neu.models.Rdata")

loadpdf <- function
### Load a data set and plot it.
(desc="",
### Data set description.
 viewer="xpdf"
### program to view output.
 ){
  if(desc!="")desc <- paste('.',desc,sep="")
  lfile <- paste("sims",desc,".models.Rdata",sep="")
  simsname <- paste("sims",desc,sep="")
  modelsname <- paste("models",desc,sep="")
  load(file=lfile)
  sims <- get(simsname)
  models <- get(modelsname)
  df <- ppp.df(sims,models)
  subt <- deduce.param.label(do.call("cbind",sims[[1]]$p[display.params]))
  print(subt)

  makepdf <- function(fun.name,arg,...){
    plot.fun <- get(fun.name)
    outf <- paste(fun.name,desc,".pdf",sep="")
    pdf(outf,h=14,w=10)
    print(plot.fun(arg,sub=subt,...))
    dev.off()
    cmd <- paste(viewer,outf,"&")
    system(cmd)
  }

  makepdf("dens.several.s",df)
  makepdf("cutoff.plot",classify.loci(df),dens=df)

  df
}

load("/home/thocking/nicholsonppp/pkg/etc/sims.neu.models.Rdata")
df.neu <- ppp.df(sims.neu,models.neu)
cl.neu <- classify.loci(df.neu)
pdf("cutoff.plot.neu.pdf",h=14,w=10,paper="a4")
cutoff.plot(cl.neu,dens=df.neu,ylim=c(-5,21)/10)
dev.off()
 
## loadpdf()
## loadpdf("few")
## loadpdf("neu",ymax=2)
foo <- classify.loci(ppp <- ppp.df(sims,models))
cutoff.plot(foo)
pdf("cutoff.pdf",paper="a4",h=0,w=0)
cutoff.plot(foo,dens=ppp)
dev.off()
system(paste("xpdf","cutoff.pdf"))
system(paste("lpr","cutoff.pdf"))


a <- mdply(data.frame(desc=c("","few","neu"),viewer="xpdf"),loadpdf)
levels(a$desc) <- c("12 populations, 1000 loci",
                    "4 populations, 1000 loci",
                    "12 populations, 19999 loci")
acl <- classify.loci(a)
cutoff.plot(acl)

pdf("roc-desc.pdf",h=6,w=8.5)
roc.loci(acl,layout=c(3,1),aspect=1)
dev.off()
system(paste("xpdf","roc-desc.pdf"))

pdf("roc-s.pdf",h=0,w=0,paper="a4")
xyplot(sensitivity~1-specificity|s,acl,
       type='l',
       groups=desc,
       panel=function(...){
         panel.abline(0,1,col="grey")
         panel.xyplot(...)
       },
       main="ROCs vary with selection strength and number of populations",
       auto.key=list(title="Simulation",space="right",lines=TRUE,points=FALSE),
       layout=c(1,5),
       aspect=1)
dev.off()