Vẽ độ lệch hồ sơ cho phù hợp GLM trong R

Question

Tôi muốn có thể vẽ sơ đồ độ lệch tiểu sử cho một ước tính tham số được trang bị sử dụng hàm glm() trong R. Sơ đồ Deviance là hàm lệch cho các giá trị khác nhau của ước tính tham số trong câu hỏi, sau khi ước tính tất cả các thông số khác. Tôi cần âm mưu sai lệch đối với một số giá trị xung quanh thông số được trang bị, để kiểm tra giả thiết hàm lệch phương trình bậc hai.

Mô hình của tôi dự đoán việc tái phạm người phạm tội. Công thức có dạng: reconv ~ [other variables] + sex, trong đó reconv là hệ số có/không nhị phân và sex là hệ số nam/nữ nhị phân. Tôi muốn vẽ độ lệch tiểu sử của tham số ước tính cho giới tính = nữ (giới tính = nam là mức tham chiếu).

Hàm glm() ước tính tham số là -0,22, với lỗi std 0,12.

[Tôi đang đặt câu hỏi này vì không có câu trả lời tôi có thể tìm thấy, nhưng tôi đã làm việc đó và muốn đăng một giải pháp để sử dụng cho người khác. Tất nhiên, sự giúp đỡ bổ sung được hoan nghênh. :-)]

Answer 1

Xem profileModel package bởi Ioannis Kosmidis. Ông đã có một bài báo trong R Journal (R News nó sẽ xuất hiện) minh họa các gói:

Ioannis Kosmidis. Gói profilemodel R: Định hướng các mục tiêu cho các mô hình với các dự báo tuyến tính. R Tin tức, 8 (2): 12-18, tháng 10 năm 2008.

PDF là here (toàn bộ bản tin).

Answer 2

Xem ?profile.glm (và example("profile.glm")) trong gói MASS - Tôi nghĩ rằng nó sẽ làm tất cả những gì bạn muốn (điều này không được nạp theo mặc định, nhưng nó được đề cập trong ?profile, mà có thể đã là nơi đầu tiên bạn nhìn ...) (Lưu ý rằng các cấu hình thường được vẽ trên một thang tỷ lệ hình vuông đã ký, để một hồ sơ bậc hai thực sự sẽ xuất hiện như một đường thẳng.)

Answer 3

Cách tôi tìm thấy để thực hiện việc này liên quan đến việc sử dụng hàm offset() (được mô tả chi tiết trong Pawitan, Y. (2001) 'In All Likelihood' p172). Các câu trả lời được đưa ra bởi @BenBolker và @GavinSimpson là tốt hơn so với điều này, trong đó họ tham chiếu các gói mà sẽ làm tất cả mọi thứ này và nhiều hơn nữa. Tôi đăng bài này bởi vì cách khác của nó làm tròn nó, và cũng có thể, vẽ đồ vật "thủ công" đôi khi rất hay cho việc học. Nó dạy tôi rất nhiều.

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