r - 在 R 中绘制轮廓似然曲线

Question

我想弄清楚如何绘制 GLM 的轮廓似然曲线
在同一图上具有 95% pLCI 的参数。我一直在尝试的例子
与下面。我得到的图不是我得到的似然曲线
期待。图的 y 轴是 tau，我想要那个轴
是可能性，以便我有一条在参数处最大的曲线
估计。我不确定在哪里可以找到这些似然值？我可能只是
误解了这背后的理论。谢谢你提供的所有帮助。

最大限度

clotting <- data.frame(
u = c(5,10,15,20,30,40,60,80,100),
lot1 = c(118,58,42,35,27,25,21,19,18),
lot2 = c(69,35,26,21,18,16,13,12,12))
glm2<-glm(lot2 ~ log(u), data=clotting, family=Gamma)
prof<-profile(glm2)
plot(prof) 

Answer 1

重新生成您的示例:

clotting <- data.frame(
                       u = c(5,10,15,20,30,40,60,80,100),
                       lot1 = c(118,58,42,35,27,25,21,19,18),
                       lot2 = c(69,35,26,21,18,16,13,12,12))
glm2 <- glm(lot2 ~ log(u), data=clotting, family=Gamma)

profile.glm函数实际上存在于 MASS包裹:

library(MASS)
prof<-profile(glm2)

为了弄清楚什么profile.glm和 plot.profile正在做，见?profile.glm和 ?plot.profile .然而，为了深入挖掘profile对象 检查 MASS:::profile.glm 的代码也可能有用和 MASS:::plot.profile ...基本上，这些告诉你的是profile返回偏差与最小偏差之差的有符号平方根，由离散参数缩放。这样做的原因是，完美二次曲线的轮廓将显示为一条直线(检测与直线的偏差比肉眼检测抛物线的偏差要容易得多)。

另一个可能有用的信息是配置文件的存储方式。基本上，它是一个数据框列表(每个参数配置一个)，除了单个数据框有点奇怪(包含一个向量分量和一个矩阵分量)。

> str(prof)
List of 2
 $ (Intercept):'data.frame':    12 obs. of  3 variables:
  ..$ tau     : num [1:12] -3.557 -2.836 -2.12 -1.409 -0.702 ...
  ..$ par.vals: num [1:12, 1:2] -0.0286 -0.0276 -0.0267 -0.0258 -0.0248 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "(Intercept)" "log(u)"
  ..$ dev     : num [1:12] 0.00622 0.00753 0.00883 0.01012 0.0114 ...
 $ log(u)     :'data.frame':    12 obs. of  2 variables:
  ..$ tau     : num [1:12] -3.516 -2.811 -2.106 -1.403 -0.701 ...
  ..$ par.vals: num [1:12, 1:2] -0.0195 -0.0204 -0.0213 -0.0222 -0.023 ...
  .. ..- attr(*, "dimnames")=List of 2

它还包含属性 summary和 original.fit您可以使用它来恢复色散和最小偏差:

disp <- attr(prof,"summary")$dispersion
mindev <- attr(prof,"original.fit")$deviance

现在反转参数 1 的转换:

dev1 <- prof[[1]]$tau^2
dev2 <- dev1*disp+mindev

阴谋:

plot(prof[[1]][,1],dev2,type="b")

(这是偏差图。您可以乘以 0.5 以获得负对数似然，或乘以 -0.5 以获得对数似然......)

编辑 : 一些更通用的函数将配置文件转换为 lattice 的有用格式/ggplot策划...

tmpf <- function(x,n) {
    data.frame(par=n,tau=x$tau,
               deviance=x$tau^2*disp+mindev,
               x$par.vals,check.names=FALSE)
}
pp <- do.call(rbind,mapply(tmpf,prof,names(prof),SIMPLIFY=FALSE))
library(reshape2)
pp2 <- melt(pp,id.var=1:3)
pp3 <- subset(pp2,par==variable,select=-variable)

现在用格子绘制它:

library(lattice)
xyplot(deviance~value|par,type="b",data=pp3,
       scales=list(x=list(relation="free")))

或者使用 ggplot2:

library(ggplot2)
ggplot(pp3,aes(value,deviance))+geom_line()+geom_point()+
    facet_wrap(~par,scale="free_x")