6

我试图弄清楚如何在同一图上绘制具有 95% pLCI 的 GLM 参数的轮廓似然曲线。我一直在尝试的例子如下。我得到的图不是我期望的似然曲线。图的 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)
4

2 回答 2

9

重新生成您的示例：

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.glmMASS:::plot.profileprofile

可能有用的另一件事是如何存储配置文件。基本上，它是一个数据帧列表（每个参数配置一个），除了单个数据帧有点奇怪（包含一个向量分量和一个矩阵分量）。

> 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")

在此处输入图像描述

于 2012-08-11T18:01:15.950 回答
1

仅供参考，为了好玩，我将上面的内容组合成一个函数，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")

在此处输入图像描述

于 2018-10-03T14:26:02.263 回答