Travis Ingram commented that it would be nice to be able to plot density maps using the function densityMap in grayscale. He gave a line of code that could be modified internally to do this; however (as of today!) this is not necessary as we can instead modify our object of class "densityMap" (now returned invisibly by the function) and replot using plot.densityMap. This is what that would look like:
> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.60’
> # this is just to simulate some data
> tree<-pbtree(n=80,scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-c(0,1)
> tree<-sim.history(tree,Q)
> x<-tree$states
> mtrees<-make.simmap(tree,x,nsim=100)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
0 1
0 -1.25418 1.25418
1 1.25418 -1.25418
(estimated using likelihood);
and (mean) root node prior probabilities
pi =
0 1
0.5 0.5
Done.
> # we don't care about this plot
> maps<-densityMap(mtrees,res=500)
sorry - this might take a while; please be patient
> # now let's change our colormap using Travis's code
> maps$cols[]<-grey(seq(1,0,length.out=length(maps$cols)))
> plot(maps,fsize=c(0.6,1),outline=TRUE,lwd=5)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.60’
> # this is just to simulate some data
> tree<-pbtree(n=80,scale=1)
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-c(0,1)
> tree<-sim.history(tree,Q)
> x<-tree$states
> mtrees<-make.simmap(tree,x,nsim=100)
make.simmap is sampling character histories conditioned on the transition matrix
Q =
0 1
0 -1.25418 1.25418
1 1.25418 -1.25418
(estimated using likelihood);
and (mean) root node prior probabilities
pi =
0 1
0.5 0.5
Done.
> # we don't care about this plot
> maps<-densityMap(mtrees,res=500)
sorry - this might take a while; please be patient
> # now let's change our colormap using Travis's code
> maps$cols[]<-grey(seq(1,0,length.out=length(maps$cols)))
> plot(maps,fsize=c(0.6,1),outline=TRUE,lwd=5)
(Click for higher res version.) Cool.
Thanks to Travis for the great suggestion & code.
this is quite cool! Not sure if you've gone over this somewhere else, but you can pretty much integrate among two (or more) colors to use any scheme you want (such as the colorbrewer ones) using the colorRampPalette function:
ReplyDeletefor example, a black-red palette similar to what plotSimmap uses:
rampfxn <- colorRampPalette(c(rgb(0,0,0), rgb(1,0,0)))
maps$cols[]<-rampfxn(length(maps$cols))
I used two colors for that ramp, but you could use any number, and the function interpolates across them. So for example, using colorbrewer:
rampfxn <- colorRampPalette(brewer.pal(11, 'RdGy'))
maps$cols[]<-rampfxn(length(maps$cols))
which can quickly get pretty crazy :)