在 R 中绘制轮廓似然曲线答案

【问题标题】：Plotting profile likelihood curves in R在 R 中绘制轮廓似然曲线
【发布时间】：2012-08-11 16:00:19
【问题描述】：

我试图弄清楚如何绘制 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)

【问题讨论】：

您从未将调用 glm 的结果存储在此处。
不知道怎么复制过来的，谢谢。

标签： r

【解决方案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")

【讨论】：

真的很有帮助，正是我想要的。非常感谢，本。
注，我认为 tau 在最近的版本中已更改为 z。

【解决方案2】：

仅供参考，为了好玩，我将上面的内容并使用purrr::imap_dfr 将其组合成一个函数，因为我找不到实现上述内容的包。

get_profile_glm <- function(aglm){
  prof <- MASS:::profile.glm(aglm)
  disp <- attr(prof,"summary")$dispersion

  purrr::imap_dfr(prof, .f = ~data.frame(par = .y,
                                deviance=.x$z^2*disp+aglm$deviance, 
                                values = as.data.frame(.x$par.vals)[[.y]],
                                stringsAsFactors = FALSE))

}

效果很好！

counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())

ggplot(get_profile_glm(aglm), aes(x = values, y = deviance)) +
  geom_point() +
  geom_line() +
  facet_wrap(~par, scale = "free_x")

【讨论】：