Wyznaczenie odchylenia profilu dla dopasowania GLM w R

Question

Chciałbym móc wykreślić odchylenie profilu dla oszacowania parametru dopasowanego za pomocą funkcji glm() w R. Odchylenie profilu jest funkcją odchylenia dla różnych wartości oszacowania parametru o której mowa, po oszacowaniu wszystkich innych parametrów. Muszę wykreślić odchylenie dla kilku wartości wokół dopasowanego parametru, aby sprawdzić założenie kwadratowej funkcji odchylenia.

Mój model przewiduje pogodzenie przestępców. Formuła ma postać: reconv ~ [other variables] + sex, gdzie reconv jest binarnym współczynnikiem tak/nie, a sex jest binarnym czynnikiem męskim/żeńskim. Chciałbym wykreślić odchylenie profilu parametru oszacowanego dla płci = kobieta (płeć = samiec jest poziomem odniesienia).

Funkcja glm() oszacowała parametr jako -0.22, z błędem std 0.12.

[pytam to pytanie, ponieważ nie było odpowiedzi udało mi się znaleźć, ale poradziłem sobie i chciał umieścić rozwiązanie być przydatne dla innych. Oczywiście dodatkowa pomoc jest mile widziana. :-)]

Answer 1

Zobacz profileModel package autorstwa Ioannis Kosmidis. Miał papier w Dzienniku R (pojawił się komunikat R) ilustrujący pakiet:

Ioannis Kosmidis. Pakiet profilemodel R: Cele profilowania dla modeli z predyktorami liniowymi. R Wiadomości, 8 (2): 12-18 października 2008.

Plik PDF jest here (cały biuletyn).

Answer 2

Zobacz ?profile.glm (i example("profile.glm")) w pakiecie MASS - myślę, że będzie robić wszystko, co chcesz (to nie jest ładowany domyślnie, ale nie są one wymienione w ?profile, co może być pierwszym miejscem wyglądasz ...) (Zwróć uwagę, że profile są zazwyczaj drukowane na podstawie skala z kwadratem kwadratowym, aby prawdziwie kwadratowy profil wyglądał jak linia prosta.)

Answer 3

Sposób znalazłem to zrobić polega na wykorzystaniu offset() funkcja (jak wyszczególniono w Pawitan, Y. (2001) według wszelkiego prawdopodobieństwa "p172). Odpowiedzi udzielone przez @BenBolker i @GavinSimpson są lepsze niż te, ponieważ odnoszą się do pakietów, które zrobią wszystko, co robi i wiele więcej. Publikuję to, ponieważ jest to w inny sposób, a także, plotkowanie rzeczy "ręcznie" jest czasami przyjemne do nauki. Wiele mnie nauczyło.

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