GLMフィットのプロファイル偏差をR

Question

にプロットするRの関数glm()を使用してフィッティングされたパラメータ推定値のプロファイル偏差をプロットすることができます。プロファイルDevianceは、パラメータ推定値のさまざまな値に対する偏差関数です他のすべてのパラメータを推定した後に、私は、2次の偏差関数の仮定をチェックするために、適合パラメータの周りのいくつかの値の偏差をプロットする必要があります。

私のモデルは犯罪者の再確認を予測しています。式は reconv ~ [other variables] + sexの形式であり、reconvはバイナリのyes/noファクタであり、sexはバイナリの男性/女性係数です。私は、性別=女性（性別=男性は参照レベル）について推定されたパラメータのプロファイル偏差をプロットしたいと思う。

glm()関数は、パラメータを-0.22と見積もり、stdエラー0.12を返しました。

[私が見つけることができる答えがなかったのでこの質問をしていますが、私はそれを試してみて、他の人に使いたい解決策を投稿したかったのです。もちろん、追加のヘルプは大歓迎です。 :-)]

Answer 1

Ioannis KosmidisによってprofileModel packageを参照してください。彼はパッケージを説明しているRジャーナル（R Newsに掲載されている）の論文を持っていた。

Ioannis Kosmidis。 profilemodel Rパッケージ：線形予測子を持つモデルの目的をプロファイリングする。 R News、8（2）：12-18、October 2008.

PDFはhere（ニュースレター全体）です。

Answer 2

MASSパッケージに?profile.glm（およびexample("profile.glm")）を参照してください - 私はそれはあなたが（これはデフォルトではロードされませんしたい全力を尽くすだろうと思うが、それは?profileに記載されている、最初の場所だった可能性がありますあなたは見ました...）（プロファイルは一般に符号付き平方根スケールでプロットされるので、真の二次プロファイルは直線として現れます）

Answer 3

これを行う方法は、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