2

I'm using the function gammamixEM from the package mixtools. How can I return the graphical output of density as in the function normalmixEM (i.e., the second plot in plot(...,which=2)) ?

Update:

Here is a reproducible example for the function gammamixEM:

x <- c(rgamma(200, shape = 0.2, scale = 14), rgamma(200, 
     shape = 32, scale = 10), rgamma(200, shape = 5, scale = 6))
out <- gammamixEM(x, lambda = c(1, 1, 1)/3, verb = TRUE)

Here is a reproducible example for the function normalmixEM:

data(faithful)
attach(faithful)
out <- normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03)
plot(out, which=2)

enter image description here

I would like to obtain this graphical output of density from the function gammamixEM.

Hack-R
  • 22,422
  • 14
  • 75
  • 131
Marine
  • 521
  • 1
  • 4
  • 23

1 Answers1

1

Here you go.

out <- normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03)

x          <- out 
whichplots <- 2
density = 2 %in% whichplots
loglik  = 1 %in% whichplots

def.par    <- par(ask=(loglik + density > 1), "mar") # only ask and mar are changed
mix.object <- x


k    <- ncol(mix.object$posterior)
x    <- sort(mix.object$x)
a    <- hist(x, plot = FALSE)
maxy <- max(max(a$density), .3989*mix.object$lambda/mix.object$sigma)

I just had to dig into the source code of plot.mixEM

So, now to do this with gammamixEM:

x <- c(rgamma(200, shape = 0.2, scale = 14), rgamma(200, 
                                                    shape = 32, scale = 10), rgamma(200, shape = 5, scale = 6))
gammamixEM.out <- gammamixEM(x, lambda = c(1, 1, 1)/3, verb = TRUE)



mix.object <- gammamixEM.out

k    <- ncol(mix.object$posterior)
x    <- sort(mix.object$x)
a    <- hist(x, plot = FALSE)
maxy <- max(max(a$density), .3989*mix.object$lambda/mix.object$sigma)

main2 <- "Density Curves"
xlab2 <- "Data" 
col2  <- 2:(k+1) 

hist(x, prob = TRUE, main = main2, xlab = xlab2, 
     ylim = c(0,maxy))


for (i in 1:k) {
  lines(x, mix.object$lambda[i] * 
          dnorm(x, 
                sd = sd(x)))
}

enter image description here

I believe it should be pretty straight forward to continue this example a bit, if you want to add the labels, smooth lines, etc. Here's the source of the plot.mixEM function.

Hack-R
  • 22,422
  • 14
  • 75
  • 131
  • Thank you very much Hack-R for your answer. How should I change the code line in `plot.mixEM`: `lines(x, mix.object$lambda[i] * dnorm(x, mean = mix.object$mu[i *arbmean + (1 - arbmean)], sd = mix.object$sigma[i *arbvar + (1 - arbvar)]), col = col2[i], lwd = lwd2)` to use the function `gammamixEM` ? Thanks a lot – Marine Aug 08 '16 at 17:25
  • @Marine Oh, ok. Now I understand what you want. I think we got lucky because the code I have now works with `gammamixEM`. I am about to update my answer right now. – Hack-R Aug 08 '16 at 17:33
  • Exactly ! I would like to plot the output of `gammamixEM` rather than `normalmixEM`. – Marine Aug 08 '16 at 17:33
  • @Marine Great. Just refresh your browser now so that you can see the new version of my answer. – Hack-R Aug 08 '16 at 17:36
  • Thank you very much Hack-R for the code. I also would like to add the fitted mixture density. I tested this answer: [link] http://stackoverflow.com/questions/17450591/fit-gamma-mixture-to-fertility-schedule-in-r. But It doesn't work in my case. – Marine Aug 08 '16 at 18:13
  • @Marine OK, I updated it with my best guess at that. To validate this or get support on obtaining this value (since `gammamixEM` doesn't have a built in plot method and lacks several of the other curve-related values from other functions like `mu` and `sigma`) I would suggest to ping CrossValidated for validation and the author's GitHub to request the feature and to see if this (or something else) is a valid workaround. – Hack-R Aug 08 '16 at 18:28