5

我希望能够绘制使用glm()R 中的函数拟合的参数估计的轮廓偏差。轮廓偏差是在估计所有其他参数之后，所讨论的参数估计的不同值的偏差函数。我需要围绕拟合参数绘制几个值的偏差，以检查二次偏差函数的假设。

我的模型是预测罪犯的再定罪。该公式的形式为： reconv ~ [other variables] + sex，其中reconv是二进制是/否因子，并且sex是二进制男性/女性因子。我想绘制为 sex=female 估计的参数的轮廓偏差（sex=male 是参考水平）。

该glm()函数将参数估计为 -0.22，标准误差为 0.12。

[我问这个问题是因为我找不到答案，但我解决了问题，并想发布一个对其他人有用的解决方案。当然，欢迎提供额外的帮助。:-)]

4

3 回答 3

6

请参阅Ioannis Kosmidis 的profileModel包。他在 R 杂志（R News 会出现）上发表了一篇论文，说明了这个包裹：

扬尼斯·科斯米迪斯。profilemodel R 包：具有线性预测变量的模型的分析目标。R 新闻，8(2):12-18，2008 年 10 月。

PDF 在这里（完整的时事通讯）。

于 2012-03-20T14:05:17.437 回答
5

请参阅包中的?profile.glm(and example("profile.glm")) MASS-- 我认为它会执行您想要的所有操作（默认情况下未加载，但在 中提到?profile，这可能是您查看的第一个地方...）（请注意，配置文件是通常绘制在有符号平方根刻度上，因此真正的二次轮廓将显示为一条直线。）

于 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 回答