6

同じプロットに95％のpLCIを持つGLMパラメーターのプロファイル尤度曲線をプロットする方法を理解しようとしています。私が試した例は以下のとおりです。私が得ているプロットは、私が期待していた尤度曲線ではありません。プロットのy軸はタウであり、パラメーター推定値で最大になる曲線が得られるように、その軸を尤度にしたいと思います。それらの尤度値がどこにあるかわかりませんか？私はこの背後にある理論を誤解しているだけかもしれません。あなたが与えることができるどんな助けにも感謝します。

マックス

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.glmありMASS:::plot.profileます...基本的に、これらが示すことは、逸脱度と最小逸脱度の差の符号付き平方根profileを返すことです。分散パラメータ。これが行われる理由は、完全な2次プロファイルのプロファイルが直線として表示されるようにするためです（放物線からの逸脱を目で検出するよりも、直線からの逸脱を検出する方がはるかに簡単です）。

知っておくと便利なもう1つのことは、プロファイルの保存方法です。基本的に、これはデータフレームのリストです（プロファイルされたパラメーターごとに1つ）。ただし、個々のデータフレームは少し奇妙です（1つのベクトルコンポーネントと1つのマトリックスコンポーネントを含みます）。

> 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参考までに、上記を実装するパッケージが見つからなかったので、楽しみのために、上記を1つの関数にまとめました。

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 に答える