library(qtl)
library(qtlMprob)

set.seed(123235823)
data(listeria)
listeria$pheno$surv264 <- ifelse(listeria$pheno$T264==264,1,0)
naff <- sum(listeria$pheno$surv264==1, na.rm=T)
nunaff <- sum(listeria$pheno$surv264==0, na.rm=T)
mymap <- pull.map(listeria)

#generate data sets
listeria.aff <- subset(listeria, ind=listeria$pheno$surv264==1)
both.index <- c(which(listeria$pheno$surv264==1), which(listeria$pheno$surv264==0)[sample(1:nunaff,naff)]) #ind for both
listeria.both <- subset(listeria, ind=both.index)
listeria.both.missing <- listeria
for(chr in 1:20)  listeria.both.missing$geno[[chr]]$data[-both.index,] <- NA

listeria <- calc.genoprob(listeria, step=1)
listeria.aff <- calc.genoprob(listeria.aff, step=1)
listeria.both <- calc.genoprob(listeria.both, step=1)
listeria.both.missing <- calc.genoprob(listeria.both.missing,step=1)

#genome scans
listeria.scan <- scanone(listeria, pheno=3, chr=1:19, model="binary") #full genotypes
maxlod.full <- max(listeria.scan$lod)

listeria.aff.scan <- scanone.reverse(listeria.aff, pheno=3, chr=1:19, phi.null="seg", step=1) #affected only
maxlod.aff <- max(listeria.aff.scan$lod)

listeria.both.scan <- scanone.reverse(listeria.both, pheno=3, chr=1:19, step=1, phi.null="em") #both, maximize denom
maxlod.both <- max(listeria.both.scan$lod)
listeria.both.scan2 <- scanone.reverse(listeria.both, pheno=3, chr=1:19, step=1, phi.null="seg", constrain=F) #both, set denom
maxlod.both2 <- max(listeria.both.scan2$lod)

listeria.both.rscan <- scanone.full(listeria.both.missing, pheno=3, chr=1:19, step=1, suppress.warnings=T) #both, restricted
maxlod.bothr <- max(listeria.both.rscan$lod)