CODEWITHSUNDEEP
rplotgamma-distribution
viridius
viridius

Reputation: 497

Add fitted (gamma) distribution density curve to histogram plot from library(MASS) fitdistr

I would like to plot the gamma density function derived from a set of observations over the histogram of the observed data. I am able to generate the histogram and parameter estimates for the gamma fit. This is being done for multiple subsets of data from a master data set. How can I plot the gamma density function on each of the histograms created in this loop?

I currently have:

library(MASS)

species <- c('acesac', 'acesac', 'acesac', 'acesac', 'acesac', 'acesac',
 'polbif', 'polbif', 'polbif', 'polbif', 'polbif', 'polbif')
tmean <- c(2,3,5,6,6,7,5,6,6,6,8,9) 
Data <- data.frame(species, tmean) 

for (i in unique(Data$species)){
  subdata <- subset(Data, species ==i)
  hist(subdata$tmean, main = i)
  dist <- fitdistr(subdata$tmean, "gamma")
}

I'm thinking that I should use lines(), however, not sure how to specify this?

Upvotes: 1

Views: 2631

Answers (1)

JasonAizkalns
JasonAizkalns

Reputation: 20473

I would add library(MASS) to your example. You may want to try doing something with curve and using add = TRUE. Another option would be to use the library(fitdistrplus) as it can plot the output of dist directly; however, I could not find a way (quickly) to change the plot titles.

library(MASS)
species <- c('acesac', 'acesac', 'acesac', 'acesac', 'acesac', 'acesac',
             'polbif', 'polbif', 'polbif', 'polbif', 'polbif', 'polbif')
tmean <- c(2,3,5,6,6,7,5,6,6,6,8,9) 
Data <- data.frame(species, tmean) 

for (i in unique(Data$species)) {
  subdata <- subset(Data, species ==i)
  dist <- fitdistr(subdata$tmean, "gamma")
  hist(subdata$tmean, main = i)
  curve(dgamma(x, shape = dist$estimate[1], rate = dist$estimate[2]), 
        add = TRUE,
        col = "red")
}

Per my comment about library(fitdistrplus) see the output from:

library(fitdistrplus)

for (i in unique(Data$species)) {
  subdata <- subset(Data, species ==i)
  dist <- fitdistrplus::fitdist(subdata$tmean, "gamma")
  plot(dist)
}

Notice that you get additional graphs (Q-Q plot, empirical and theoretical CFS, and P-P plot) and the density lines are plotted for "free". However, you lose the ability to add main = i. I'm sure someone smarter than me can figure out a quick way to add titles or modify the fitdistrplus::plot.fitdist method - might be worth asking a separate question.

Upvotes: 1

Related Questions