Take a glm and creates a deviance profile, as MASS::profile only returns a transformed version of the deviance.

glm_profile.R

library(ggplot2)

#A glm from the glm helpfile
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())

#The old way...
prof <- MASS:::profile.glm(glm.D93)
MASS:::plot.profile(prof)


#some info from https://stackoverflow.com/questions/11916139/plotting-profile-likelihood-curves-in-r
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))
  
}


#make and plot a profile
ggplot(get_profile_glm(glm.D93), aes(x = values, y = deviance)) +
  geom_point() +
  geom_line() +
  facet_wrap(~par, scale = "free_x")


#Or with another package
library(profileModel)
prof_mod <- profileModel(glm.D93, quantile = qchisq(0.95, 1) ,
             objective = "ordinaryDeviance")

#nice
plot(prof_mod)
profConfint(prof_mod)
pairs(prof_mod)

#ggplot2?