5

R の関数を使用してフィッティングされたパラメーター推定値のプロファイル逸脱度をプロットできるようにしたいと考えていますglm()。プロファイル逸脱度は、他のすべてのパラメーターを推定した後の、問題のパラメーター推定値のさまざまな値の逸脱度関数です。二次逸脱関数の仮定を確認するために、適合パラメータの周りのいくつかの値の逸脱をプロットする必要があります。

私のモデルは、犯罪者の再有罪判決を予測しています。式の形式は次の reconv ~ [other variables] + sexとおりreconvですsex。性別=女性 (性別=男性は参考レベル) について推定されたパラメーターのプロファイル逸脱度をプロットしたいと思います。

関数はglm()パラメーターを -0.22 と推定し、標準誤差は 0.12 でした。

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

4

3 に答える 3

6

Ioannis Kosmidisによる profileModelパッケージを参照してください。彼は R Journal (R News と表示される) にパッケージを説明する論文を持っていました。

イオアニス・コスミディス profilemodel R パッケージ: 線形予測子を持つモデルのプロファイリング目標。R ニュース、8(2):12-18、2008 年 10 月。

PDFはこちら（ニュースレター全体）。

于 2012-03-20T14:05:17.437 に答える
5

?profile.glmパッケージ内の(およびexample("profile.glm")) を参照してくださいMASS-- 必要なことはすべて実行されると思います (これはデフォルトではロードされ?profileませんが、最初に見た場所である可能性がある に記載されています ...) (プロファイルは一般に、正の 2 次プロファイルが直線として表示されるように、符号付き平方根スケールでプロットされます。)

于 2012-03-20T13:58:56.300 に答える
0

これを行うために私が見つけた方法には、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
于 2012-03-23T16:10:58.793 に答える