Restore your example:
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)
The profile.glm function actually lives in the MASS package:
library(MASS) prof<-profile(glm2)
To find out what profile.glm and plot.profile , see ?plot.profile and ?plot.profile . However, in order to dig up the profile object, it may also be useful to examine the code MASS:::profile.glm and MASS:::plot.profile ... basically, what they tell you is that profile returns the square root of the difference between the deviations and minimum deviation, scalable by the variance parameter. The reason this is done is because the profile of the perfectly quadratic profile will be displayed as a straight line (it is much easier to detect deviations from a straight line than from a parabola in the eye).
Another thing that may be useful to know is how to save a profile. Basically, this is a list of data frames (one for each parameter of the profiled), except that the individual data frames are a bit strange (contain one vector component and one matrix component).
> str(prof) List of 2 $ (Intercept):'data.frame': 12 obs. of 3 variables: ..$ tau : num [1:12] -3.557 -2.836 -2.12 -1.409 -0.702 ... ..$ par.vals: num [1:12, 1:2] -0.0286 -0.0276 -0.0267 -0.0258 -0.0248 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] "(Intercept)" "log(u)" ..$ dev : num [1:12] 0.00622 0.00753 0.00883 0.01012 0.0114 ... $ log(u) :'data.frame': 12 obs. of 2 variables: ..$ tau : num [1:12] -3.516 -2.811 -2.106 -1.403 -0.701 ... ..$ par.vals: num [1:12, 1:2] -0.0195 -0.0204 -0.0213 -0.0222 -0.023 ... .. ..- attr(*, "dimnames")=List of 2
It also contains the attributes summary and original.fit , which you can use to restore variance and minimize deviation:
disp <- attr(prof,"summary")$dispersion mindev <- attr(prof,"original.fit")$deviance
Now change the conversion for parameter 1:
dev1 <- prof[[1]]$tau^2 dev2 <- dev1*disp+mindev
Plot:
plot(prof[[1]][,1],dev2,type="b")
(This is the deviation graph. You can multiply by 0.5 to get negative log-likelihood, or -0.5 to get log-likelihood ...)
edit : some more general functions for converting a profile into a useful format for lattice / ggplot plotting ...
tmpf <- function(x,n) { data.frame(par=n,tau=x$tau, deviance=x$tau^2*disp+mindev, x$par.vals,check.names=FALSE) } pp <- do.call(rbind,mapply(tmpf,prof,names(prof),SIMPLIFY=FALSE)) library(reshape2) pp2 <- melt(pp,id.var=1:3) pp3 <- subset(pp2,par==variable,select=-variable)
Now build it with a grid:
library(lattice) xyplot(deviance~value|par,type="b",data=pp3, scales=list(x=list(relation="free")))
Or with ggplot2:
library(ggplot2) ggplot(pp3,aes(value,deviance))+geom_line()+geom_point()+ facet_wrap(~par,scale="free_x")