0

I have the following data

[1] 0.09733344 0.17540020 0.14168188 0.54093074 0.78151039 0.28068527 [7] 1.96164429 0.33743328 0.05200734 0.09103039 0.28842044 0.09240131 [13] 0.09143535 0.38142022 0.11700952

from which I did bayesian inference and made a plot with the following code

f_theta<-function(theta,Data){
(theta^length(Data) )*exp(-theta*sum(Data))}
theta<-seq(1,20,length=100)
a=b=0.001
plot(theta,dgamma(theta,a,b),type="l",col="red",
ylim=c(0,2),tck=-0.01,cex.lab=0.8,cex.axis=0.8)
lines(theta,dgamma(theta,length(Data)+a,sum(Data)+b),col="green",lty=1)
lines(theta,f_theta(theta,Data=Data),lty=1,col="blue")
legend('topright',legend=c("Prior","Post","Likelihood")
 ,col=c("red","green","blue","purple"),lty=1,bty="n",cex=0.8)

But I've seen the following graph

enter image description here

which has code

# ggplot2 examples
library(ggplot2) 

# create factors with value labels 
mtcars$gear <- factor(mtcars$gear,levels=c(3,4,5),
labels=c("3gears","4gears","5gears")) 
mtcars$am <- factor(mtcars$am,levels=c(0,1),
labels=c("Automatic","Manual")) 
mtcars$cyl <- factor(mtcars$cyl,levels=c(4,6,8),
labels=c("4cyl","6cyl","8cyl")) 

# Kernel density plots for mpg
# grouped by number of gears (indicated by color)
qplot(mpg, data=mtcars, geom="density", fill=gear, alpha=I(.5), 
main="Distribution of Gas Milage", xlab="Miles Per Gallon", 
ylab="Density")

but I'm not quite familiar with ggplot library and graphs and I would like some help in order to adapt my code and make a graph similar to last one.

1
  • See if you can narrow this down to the core of the question: you have some data (you could post the output of your analysis), you want to plot distributions, and you're not sure what bits of the example ggplot code need to be changed
    – camille
    Commented May 23, 2018 at 20:51

1 Answer 1

3

ggplot() assumes that your data are in a particular format (sometimes called "long", but the author of ggplot() dislikes that description), so let's start by putting them into that format:

Data2 = data.frame(
    theta = rep(theta, 3),
    WhichDistribution = c(rep("Prior",length(theta)), rep("Post",length(theta)), rep("Likelihood",length(theta))),
    Density = c(dgamma(theta,a,b), dgamma(theta,length(Data)+a,sum(Data)+b), f_theta(theta,Data=Data))
)

Then we can construct a ggplot() command. ggplot() needs data, aesthetics, and a geometry. Your data will be the data frame just constructed. The aesthetics refer generally to how the qualities of the data will impact the graph (what is on axes, what determines groups, etc.), and the geometry is the kind of plot (not a great wording, sorry).

ggplot(Data2, aes(x=theta, y=Density, group=WhichDistribution, color=WhichDistribution, fill=WhichDistribution))+
    # position="identity" in order to not stack the densities
    geom_area(alpha=.2, position="identity") +
    # gets rid of the title on the legend
    theme(legend.title = element_blank())+
    # make the horizontal axis label pretty
    scale_x_continuous(expression(theta))

You can change alpha to adjust transparency. If you want the horizontal axis to not go all the way to 20, change it in scale_x_continuous():

ggplot(Data2, aes(x=theta, y=Density, group=WhichDistribution, color=WhichDistribution, fill=WhichDistribution))+
    # position="identity" in order to not stack the densities
    geom_area(alpha=.2, position="identity") +
    # gets rid of the title on the legend
    theme(legend.title = element_blank())+
    # make the horizontal axis label pretty
    scale_x_continuous(expression(theta), limits=c(0,7))

enter image description here

qplot() is a quick plotting function that seems to mostly get in the way for people trying to learn the ggplot() language, so you might want to avoid it.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.