Plotting profile likelihood curves in R

user_1592443 0
Viewed 4333

I am trying to figure out how to plot the profile likelihood curve of a GLM parameter with 95% pLCI's on the same plot. The example I have been trying with is below. The plots I am getting are not the likelihood curves that I was expecting. The y-axis of the plots is tau and I would like that axis to be the likelihood so that I have a curve that maxes at the parameter estimate. I am not sure where I find those likelihood values? I may just be misinterpreting the theory behind this. Thanks for any help you can give.

Max

clotting <- data.frame(
u = c(5,10,15,20,30,40,60,80,100),
lot1 = c(118,58,42,35,27,25,21,19,18),
lot2 = c(69,35,26,21,18,16,13,12,12))
glm2<-glm(lot2 ~ log(u), data=clotting, family=Gamma)
prof<-profile(glm2)
plot(prof)
r
2 Answers

FYI, for fun, I took the above and whipped it together into a single function using purrr::imap_dfr as I couldn't find a package that implements the above.

get_profile_glm <- function(aglm){
  prof <- MASS:::profile.glm(aglm)
  disp <- attr(prof,"summary")$dispersion

  purrr::imap_dfr(prof, .f = ~data.frame(par = .y,
                                deviance=.x$z^2*disp+aglm$deviance, 
                                values = as.data.frame(.x$par.vals)[[.y]],
                                stringsAsFactors = FALSE))

}

Works great!

counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())

ggplot(get_profile_glm(aglm), aes(x = values, y = deviance)) +
  geom_point() +
  geom_line() +
  facet_wrap(~par, scale = "free_x")

enter image description here

user_1903520
Related