Geek Answers Handbook

Draw a Deviance profile to match GLM in R

I would like to be able to build a deviation profile for estimating a parameter set using a function glm()in R. Deviance profile is a deviation function for different values ​​of the considered parameter estimation, after evaluating all other parameters, I need to build a deviation for several values ​​around the set parameter so that check the assumption of a quadratic deviation function.

My model predicts the reorganization of criminals. The formula has the form:, reconv ~ [other variables] + sexwhere reconvis the binary yes / no coefficient, and sexis the binary male / female factor. I would like to build the deviation profile of the parameter estimated for sex = female (sex = male is the reference level).

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

[I ask this question because there was no answer that I could find, but I worked through it and wanted to publish a solution that might be useful to others. Of course, additional help is appreciated. :-)]

+5
source share
3 answers

Under profileModel by Ioannis Kosmidis. He had an article in the R Journal (R News that appeared) illustrating the package:

Ioannis Cosmidis. Profilemodel R package: Profiling tasks for models with linear predictors. R News, 8 (2): October 12-18, 2008

PDF here (full newsletter).

+6
source

. ?profile.glm ( example("profile.glm")) MASS - , , ( , ?profile, , , , ...) ( , , .)

+5

, , offset() ( Pawitan, Y. (2001) "In All Likelihood" p172). @BenBolker @GavinSimpson , , , , , , . , , , "", . .

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
0

More articles:


All Articles