1

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?

JasonAizkalns
  • 20,243
  • 8
  • 57
  • 116
viridius
  • 477
  • 5
  • 17
  • What is `Data`? Can you provide a reproducible example? – davechilders Oct 09 '15 at 18:50
  • My apologies for not being able to provide a reproducible example. `Data` is the full dataset I am working with. The vector `subdata$tmean` is a set of temperature values in degrees C. – viridius Oct 09 '15 at 18:54
  • See [here](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) about supplying a reproducible example. You can always use `dput()`. – davechilders Oct 09 '15 at 18:55
  • [There are many ways to provide a reproducible example](http://stackoverflow.com/q/5963269/903061). Give it a go. It **shouldn't** be your full data, it should be a small, representation that illustrates the problem, based on a subset of real data or on simulated data. – Gregor Thomas Oct 09 '15 at 18:55
  • Ok, I'll give it a try... – viridius Oct 09 '15 at 18:57
  • 1
    I added an example for `Data` – viridius Oct 09 '15 at 19:09

1 Answers1

1

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.

JasonAizkalns
  • 20,243
  • 8
  • 57
  • 116
  • I added this line for `curve` (MASS was already loaded) and receive no errors, but it does not show the curve? – viridius Oct 09 '15 at 19:26
  • @user3830407 I meant adding `library(MASS)` to your original question so that it is fully reproducible. Are you sure the curve is not plotting? Try copying and pasting the exact code above. – JasonAizkalns Oct 09 '15 at 19:30
  • The curve plots fine using the example dataset, however, it does not plot when using the actual `Data` (I'll have to check my data). ...I was able to get the very nice set of plots using fitdistrplus! Many thanks for both solutions. – viridius Oct 09 '15 at 19:35
  • @user3830407 - perhaps the actual data does not follow a gamma distribution. best of luck and feel free to ask more reproducible questions. – JasonAizkalns Oct 09 '15 at 22:50