CODEWITHSUNDEEP
rplotglm
MatW
MatW

Reputation: 315

Plot the Profile Deviance for a GLM fit in R

I would like to be able to plot the profile deviance for a parameter estimate fitted using the function glm() in R. The profile Deviance is the deviance function for different values of the parameter estimate in question, after estimating all other parameters. I need to plot the deviance for several values around the fitted parameter, to check the assumption of quadratic deviance function.

My model is predicting reconviction of offenders. The formula is of the form: reconv ~ [other variables] + sex, where reconv is a binary yes/no factor, and sex is binary male/female factor. I'd like to plot the profile deviance of the parameter estimated for sex=female (sex=male is the reference level).

The glm() function estimated the parameter as -0.22, with std error 0.12.

[I'm asking this question because there was no answer I could find, but I worked it out, and wanted to post a solution to be of use to others. Of course, additional help is welcome. :-)]

Upvotes: 5

Views: 1480

Answers (3)

MatW
MatW

Reputation: 315

The way I found to do this involves using the offset() function (as detailed in Pawitan, Y. (2001) 'In All Likelihood' p172). The answers given by @BenBolker and @GavinSimpson are better than this, in that they referenced packages which will do everything this does and a lot more. I'm posting this because its another way round it, and also, plotting things "manually" is sometimes nice for learning. It taught me a lot.

sexi <- as.numeric(data.frame$sex)-1      #recode a factor as 0/1 numeric

beta <- numeric(60)              #Set up vector to Store the betas
deviance <- numeric(60)          #Set up vector to Store the deviances

for (i in 1:60){

  beta[i] <- 0.5 - (0.01*i)  
  #A vector of values either side of the fitted MLE (in this case -0.22)

  mod <- update(model,
                   .~. - sex             #Get rid of the fitted variable
                   + offset(   I(sexi*beta[i])   )   #Replace with offset term.
                )
  deviance[i] <- mod$deviance                        #Store i'th deviance
}

best <- which.min(deviance)                   
#Find the index of best deviance. Should be the fitted value from the model

deviance0 <- deviance - deviance[best]         
#Scale deviance to zero by subtracting best deviance

betahat <- beta[best]    #Store best beta. Should be the fitted value.
stderror <- 0.12187      #Store the std error of sex, found in summary(model)

quadratic <- ((beta-betahat)^2)*(1/(stderror^2))  
#Quadratic reference function to check quadratic assumption against

x11()                                    
plot(beta,deviance0,type="l",xlab="Beta(sex)",ylim=c(0,4))    
lines(beta,quadratic,lty=2,col=3)           #Add quadratic reference line
abline(3.84,0,lty=3)                #Add line at Deviance = 3.84

Upvotes: 0

Gavin Simpson
Gavin Simpson

Reputation: 174813

See the profileModel package by Ioannis Kosmidis. He had a paper in the R Journal (R News it would appear) illustrating the package:

Ioannis Kosmidis. The profilemodel R package: Profiling objectives for models with linear predictors. R News, 8(2):12-18, October 2008.

The PDF is here (entire newsletter).

Upvotes: 6

Ben Bolker
Ben Bolker

Reputation: 226192

See ?profile.glm (and example("profile.glm")) in the MASS package -- I think it will do everything you want (this is not loaded by default, but it is mentioned in ?profile, which might have been the first place you looked ...) (Note that the profiles are generally plotted on a signed-square-root scale, so that a truly quadratic profile will appear as a straight line.)

Upvotes: 5

Related Questions