#' Plot the results of a mixture model detection function fit.
#'
#' Plots the detection function (or pdf) of a fitted mixture model detection
#' function, optionally overlayed on a histogram of the observed data.
#'
#' @param x a \code{\link{ds.mixture}} object.
#' @param style If set to "comp", composite plots of the detection function will
#' be shown for the detection function (the averaged detection function in
#' the covariate case).
#' @param main (A vector of) title(s) for the plot(s). By default these are set
#' by the function (and are fairly ugly but descriptive).
#' @param breaks Breaks to be used for the histogram. This can be a vector of
#' numbers or any of the permissable options used in \code{\link{hist}}.
#' Defaults to "Sturges".
#' @param ylim Used to manually set the y limit of the plot. Defaults to NULL.
#' @param xlim Used to manually set the x limit of the plot. Defaults to NULL.
#' @param pdf Should the pdf be plotted rather than detection function be
#' plotted? Only really useful with point transect data. Defaults to FALSE.
#' @param plot.formula Formula of covariates to be plotted. Defaults to
#' \code{NULL}, which plots all covariates. No effect with non-covariate
#' models.
#' @param hide.hist Should the histogram be hidden, leaving only the detection
#' function (or pdf) to be plotted? Defaults to FALSE. Can be a vector the same length as the number of plots.
#' @param nomf Should the mfrow value be altered? Useful when creating custom
#' plots for publication. Defaults to FALSE (yes, change the mfrow value).
#' @param x.axis Set the x axis labels. Again, useful for publication plots.
#' Defaults to NULL, which uses the default R values.
#' @param xlab Label for the x axis.
#' @param ylab Label for the y axis.
#' @param ... not used at the moment
#' @return a plot!
#'
#' @section Details:
#' For covariate models, all the levels of factor variables are plotted or the
#' 25, 50 and 75th percentiles of continuous variables are plotted averaged over
#' the values of the other covariates.
#'
#' @author David L, Miller
#' @examples
#' library(mmds)
#' set.seed(0)
#' ## simulate some line transect data from a 2 point mixture
#' sim.dat<-sim.mix(c(-0.223,-1.897,inv.reparam.pi(0.3)),2,100,1)
#' ## fit the model
#' fit.sim.dat<-fitmix(sim.dat,1,2)
#' ## plot
#' plot(fit.sim.dat)
#'
#' @export
#' @S3method plot ds.mixture
#' @method plot ds.mixture
plot.ds.mixture<-function(x,style="",main="",breaks="Sturges",ylim=NULL,xlim=NULL,pdf=FALSE,plot.formula=NULL,hide.hist=FALSE,nomf=FALSE,x.axis=NULL,xlab="Distance",ylab=NULL,...){
fit.object<-x
# todo:
# * level legend
width<-fit.object$width
data<-fit.object$data
pars<-fit.object$pars
mix.terms<-fit.object$mix.terms
zdim<-fit.object$zdim
z<-fit.object$z
pt<-fit.object$pt
# plot resolution
plot.res<-1000
# set up the parameters, shortcuts
swpars<-switchpars(pars,mix.terms,z=z,zdim=zdim)
gp<-getpars(swpars$fpar,mix.terms,swpars$zdim,swpars$z)
sigma<-gp$key.scale
pis<-gp$mix.prop
# use the correct function for effective strip width
intfcn<-integrate.hn
if(pt)
intfcn<-integrate.hn.pt
# don't plot the histogram if asked
#if(hide.hist){
# hist.col<-"white"
#}else{
# hist.col<-"black"
#}
hist.col <- rep("white",length(hide.hist))
hist.col[!hide.hist] <- "black"
if(is.null(z)){
##### no covariates
# calculate mu
mu<-c()
for(j in 1:mix.terms)
mu<-c(mu,intfcn(sigma[j],width))
# create a sequence to evaluate the detection function, and evaluate
x.seq<-seq(0,width,len=plot.res)
plotvals<-matrix(0,mix.terms,plot.res)
for(j in 1:mix.terms)
plotvals[j,]<-keyfct.hn(x.seq,sigma[j])
# what is the value of the detection function at 0 distance?
p.at.zero<-c()
for(j in 1:mix.terms)
p.at.zero<-c(p.at.zero,pis[j]*keyfct.hn(0,sigma[j]))
p.at.zero<-sum(p.at.zero)
# make the histogram object
a<-hist(data$distance,plot=FALSE,breaks=breaks)
if(is.null(ylab)){
ylabel<-"Probability of detection"
}else{
ylabel<-ylab
}
# point transect data
if(pt){
# do pdf plotting
if(pdf){
plotvals<-matrix(0,mix.terms,plot.res)
for(j in 1:mix.terms)
plotvals[j,]<-2*pi*pis[j]*x.seq*keyfct.hn(x.seq,sigma[j])/
intfcn(sigma[j],width)
a$density<-a$density*(1/sum(a$density*diff(a$breaks)))
# this isn't really p at zero, it's the max on the y axis
p.at.zero<-max(plotvals,a$density)
if(is.null(ylab)){
ylabel<-"PDF of detected distances"
}else{
ylabel<-ylab
}
# not pdf plotting
}else{
# rescale hist by distance
rsf<-a$mids
musum<-sum(mu*pis)/(2*pi)
a$density<-a$density*(musum/sum(a$density*diff(a$breaks)))
a$density<-a$density/rsf
}
}else{
musum<-sum(mu*pis)
hist.area<-sum(a$density*diff(a$breaks))
a$density<-a$density*(musum/hist.area)
}
par(las=1)
# what is the upper limit of the plot?
if(hide.hist){
ylim.upper <- p.at.zero
}else{
ylim.upper <- max(a$density,p.at.zero)
}
# actually do the plotting
plot(a,freq=FALSE,ylab=ylabel,axes=FALSE,
ylim=c(0,ylim.upper),xlim=c(0,width),
xlab=xlab,border=hist.col,
main=main)
if(!is.null(x.axis)){
axis(1,x.axis)
}else{
axis(1)
}
axis(2,at=c(0,p.at.zero/2,p.at.zero),labels=c(0,0.5,1))
if(mix.terms>1){
lines(x.seq,pis%*%plotvals)
}else{
lines(x.seq,plotvals)
}
if(style == "comp" & mix.terms>1){
for(i in 1:mix.terms){
lines(x.seq,plotvals[i,],lty=2)
}
}
box()
}else{
##### covariate code
if(is.null(ylab)){
ylabel<-"Probability of detection"
}else{
ylabel<-ylab
}
# calculate mu per observation
mu<-apply(sigma,1,intfcn,width)
mus<-mu%*%matrix(pis,length(pis),1)
##### evaluate the detection function
# create a sequence of xs to evaluate the detection function at
x.seq<-seq(0,width,len=1000)
# quick detection function function
eval.detfct<-function(x,sigma,pis){
matrix(pis,1,length(pis))%*%apply(as.matrix(x),1,keyfct.hn,key.scale=sigma)
}
# evaluate the detection function at x.seq for each sigma
# 1000 x # samples matrix results
plotvals.mat<-apply(sigma,2,eval.detfct,pis=pis,x=x.seq)
if(pt & pdf){
plotvals.mat<-2*pi*x.seq*plotvals.mat/as.vector(mus)
}
# add in a weighting
pas<-mus/width
ws<-(1/pas)/sum(1/pas)
plotvals.mat.w<-t(t(plotvals.mat)*c(ws))
# average according to the weighting
plotvals<-rowSums(plotvals.mat.w)
##### what is the value of the detection function at 0 distance?
# evaluate the detection function at 0 distance
paz.mat<-matrix(pis,1,length(pis))%*%apply(sigma,2,keyfct.hn,distance=0)
if(pt & pdf){
paz.mat<-paz.mat/as.vector(mus)
}
# put in the weights
paz.mat<-t(t(paz.mat)*c(ws))
# average
p.at.zero<-rowSums(paz.mat)
##### calculate the overall mu (ie. area under the detection function)
mu<-sum(c(mus)*c(ws))
##### make the histogram object
a<-hist(data$distance,plot=FALSE,breaks=breaks)
if(pt){
# do pdf plotting
if(pdf){
a$density<-a$density*(1/sum(a$density*diff(a$breaks)))
# this isn't really p at zero, it's the max on the y axis
p.at.zero<-max(plotvals,a$density)
if(is.null(ylab)){
ylabel<-"PDF of detected distances"
}else{
ylabel<-ylab
}
style<-""
# not pdf plotting
}else{
# rescale hist by distance
rsf<-a$mids
musum<-sum(mu*pis)/(2*pi)
a$density<-a$density*(musum/sum(a$density*diff(a$breaks)))
a$density<-a$density/rsf
}
}else{
# re-weight the heights to the area under histogram == area under curve
a$density<-a$density*(mu/sum(a$density*diff(a$breaks)))
}
# now store plotvals into a matrix
# the covariates are cbinded to this...
pv<-matrix(plotvals,length(plotvals),1)
colnames(pv)[ncol(pv)]<-"Average detection function"
plot.names<-"Average detection function"
# store the plot sequence
# vector of number of "levels" for each covar so
# the graphics line up...
plot.seq<-c(1)
if(style=="comp"){
plotvals.comp<-c()
for(j in 1:mix.terms){
curr <- apply(matrix(sigma[j,],1,length(sigma[j,])),
2,keyfct.hn,distance=x.seq)
# add in a weighting
pas<-mus/width
ws<-(1/pas)/sum(1/pas)
plotvals.mat.w<-t(t(curr)*c(ws))
# average according to the weighting
plotvals.comp<-rbind(plotvals.comp,rowSums(plotvals.mat.w))
}
plot.seq[1]<-plot.seq[1]+mix.terms
pv<-cbind(pv,t(plotvals.comp))
}
##### plot covariate levels
# default to all
if(is.null(plot.formula))
plot.formula<-fit.object$model.formula
# do nothing if the user asks
if((as.character(plot.formula) != "none") & !(pt&pdf)){
# make sure it's a formula
plot.formula<-as.formula(plot.formula)
# work out levels etc...
for(curr.term in attr(terms.formula(plot.formula),"term.labels")){
# if we have a factor covariate
if(grepl("^as.factor",curr.term)){
# # indicator for column that we want
# ind <- colnames(data)==gsub("^as.factor\\((\\w+)\\)$","\\1",
# curr.term)
# # pull that column out
# this.col<-data[,ind]
# # find its levels
# this.col.levels<-levels(as.factor(
# data[,gsub("^as.factor\\((\\w+)\\)$","\\1",
# curr.term)]))
#
# for(tl in this.col.levels){
# ind0<-this.col==tl
#
# fac0<-sum(a$density*diff(a$breaks))/mean(mus[ind0])
#
# plotvals.mat0<-plotvals.mat[,ind0]
# pas0<-mus[ind0]/width
# ws0<-(1/pas0)/sum(1/pas0)
#
# plotvals.mat0<-t(t(plotvals.mat0)*c(ws0))
# plotvals0<-rowSums(plotvals.mat0)*fac0
# pv<-cbind(pv,plotvals0/plotvals0[1])
# colnames(pv)[ncol(pv)]<-paste(curr.term,"==",tl,sep="")
#
# }
# plot.names<-c(plot.names,
# paste("Levels of ",curr.term,sep=""))
# continuous covariates
# estimate the 0.25, 0.5 and 0.75 quantiles
this.col.levels<-colnames(z[[1]])[grepl(curr.term,colnames(z[[1]]),fixed=TRUE)]
}else{
this.col.levels<-quantile(z[[1]][,curr.term],c(0.25,0.5,0.75),names=FALSE)
}
# do the zero level for factors
if(grepl("^as.factor",curr.term)){
z.copy <- z
z.copy[[1]][,grepl(curr.term,colnames(z[[1]]),fixed=TRUE)] <- 0
swpars.copy<-switchpars(pars,mix.terms,z=z.copy,zdim=zdim)
gp.copy<-getpars(swpars.copy$fpar,mix.terms,swpars.copy$zdim,swpars.copy$z)
sigma.copy<-gp.copy$key.scale
plotvals.mat<-apply(sigma.copy,2,eval.detfct,pis=pis,x=x.seq)
plotvals0<-rowSums(plotvals.mat)
pv<-cbind(pv,plotvals0/plotvals0[1])
colnames(pv)[ncol(pv)]<-paste(curr.term,"==",0,sep="")
}
for(tl in this.col.levels){
z.copy <- z
if(grepl("^as.factor",curr.term)){
z.copy[[1]][,tl] <- 1
z.copy[[1]][,colnames(z[[1]])[grepl(curr.term,colnames(z[[1]]),
fixed=TRUE)]!=tl] <- 0
}else{
z.copy[[1]][,curr.term] <- rep(tl,nrow(z.copy[[1]]))
}
swpars.copy<-switchpars(pars,mix.terms,z=z.copy,zdim=zdim)
gp.copy<-getpars(swpars.copy$fpar,mix.terms,swpars.copy$zdim,swpars.copy$z)
sigma.copy<-gp.copy$key.scale
plotvals.mat<-apply(sigma.copy,2,eval.detfct,pis=pis,x=x.seq)
plotvals0<-rowSums(plotvals.mat)
pv<-cbind(pv,plotvals0/plotvals0[1])
if(grepl("^as.factor",curr.term)){
colnames(pv)[ncol(pv)]<-paste(curr.term,"==",
gsub("^as.factor\\(\\w+\\)(.*)$","\\1",curr.term),sep="")
}else{
colnames(pv)[ncol(pv)]<-paste(curr.term," ",tl," quantile",sep="")
}
}
if(grepl("^as.factor",curr.term)){
plot.names<-c(plot.names,paste("Levels of ",curr.term,sep=""))
plot.seq<-c(plot.seq,length(this.col.levels)+1)
}else{
plot.names<-c(plot.names,
paste(curr.term," 0.25/0.5/0.75 quantiles",sep=""))
plot.seq<-c(plot.seq,length(this.col.levels))
}
# ind<-colnames(data)==curr.term
# # pull that column out
# this.col<-data[,ind]
# # find its "levels" - the 0.25, 0.5, 0.75 quantiles
# this.col.levels<-quantile(this.col,c(0.25,0.5,0.75),names=FALSE)
#
# for(tl in this.col.levels){
# ind0<-this.col<=tl
#
# fac0<-sum(a$density*diff(a$breaks))/mean(mus[ind0])
# plotvals.mat0<-plotvals.mat[,ind0]
# pas0<-mus[ind0]/width
# ws0<-(1/pas0)/sum(1/pas0)
#
# plotvals.mat0<-t(t(plotvals.mat0)*c(ws0))
# plotvals0<-rowSums(plotvals.mat0)*fac0
# pv<-cbind(pv,plotvals0/plotvals0[1])
# colnames(pv)[ncol(pv)]<-paste(curr.term," ",tl," quantile",sep="")
# }
# plot.names<-c(plot.names,
# paste(curr.term," 0.25/0.5/0.75 quantiles",sep=""))
# }
# plot.seq<-c(plot.seq,length(this.col.levels))
}
# how big does mfrow have to be?
if(length(plot.seq)>4){
mfrows<-c(ceiling(length(plot.seq)/4),4)
}else{
mfrows<-c(1,length(plot.seq))
}
}else{
mfrows<-c(1,1)
}
# the user can specify the names if they want...
plot.names[1:length(main)]<-main
# save before changing!
if(!nomf){
op<-par()
par(mfrow=mfrows)
}
par(las=1)
k<-1
# histogram colours
hist.col <- rep(hist.col,length.out=length(plot.seq))
for(i in 1:length(plot.seq)){
# what is the upper limit of the plot?
if(hist.col[i]=="white"){
ylim.upper <- p.at.zero
}else{
ylim.upper <- max(a$density,p.at.zero)
}
# dummy plot to set things up
plot(x=c(0,width),
y=c(0,max(a$density,p.at.zero)),
type="n",
ylab=ylabel,axes=FALSE,
ylim=c(0,ylim.upper),
xlim=c(0,width),
xlab=xlab,
main=plot.names[i],las=1)
# plot the histogram
plot(a,freq=FALSE,add=TRUE,border=hist.col[i])
if(!is.null(x.axis)){
axis(1,x.axis)
}else{
axis(1)
}
# get the labels right
if(pt){
axis(2,at=c(0,p.at.zero/2,p.at.zero),
labels=round(c(0,p.at.zero/2,p.at.zero),3))
}else{
axis(2,at=c(0,p.at.zero/2,p.at.zero),labels=c(0,0.5,1))
}
box()
if(i==1){
line.style<-c(1,rep(2,mix.terms))
line.colours<-rep(1,plot.seq[i])
}else{
line.style<-rep(1,plot.seq[i])
line.colours<-grey(seq(0,0.5,len=plot.seq[i]))
}
for(j in 1:plot.seq[i]){
lines(x.seq,pv[,k],col=line.colours[j],lty=line.style[j])
k<-k+1
}
}
# restore pars
if(!nomf){
par(op)
}
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.