如何繪製物流GLM預測的R值和置信區間

我有一個存在/不存在響應變量的二項式glm並用9個電平這樣的因素變量：如何繪製物流GLM預測的R值和置信區間

data$y<-factor(data$y,levels=c(0,1),labels=c("absent","present")) 
table(data$y,data$site_name) 

      Andulay Antulang Basak Dauin Poblacion District 1 Guinsuan Kookoo's Nest Lutoban Pier Lutoban South Malatapay Pier 
    absent  4  4  1       0  3    1   5    5    2 
    present  2  2  5       6  3    5   1    1    4 

model <- glm(y~site_name,data=data,binomial)

只是跳過模型推斷和爲了簡潔起見，我如何繪製每個站點以可信區間獲取「呈現」在箱形圖中的概率？我想要的是Plot predicted probabilities and confidence intervals in R中顯示的內容，但我想用boxplot顯示它，因爲我的迴歸變量site_name是一個有9個級別的因子，而不是連續變量。

我想我可以按照以下方法計算必要的值（但不是100％肯定的正確性）：

功能的模型係數轉換回成功的概率：預測

calc_val <- function(x){return(round(1/(1+1/(exp(x))),3))}

prob <- tapply(predict(model,type="response"),data$site_name,function(x){round(mean(x),3)}) 
means <- as.data.frame(prob)

75％和95％置信區間爲預測概率：基於所述模型概率

ci <- cbind(confint(model,level=0.9),confint(model,level=0.5)) 
rownames(ci) <- gsub("site_name","",rownames(ci)) 
ci <- t(apply(ci,1,calc_val))

加入它一起在一個表

ci<-cbind(means,ci) 
ci 
          prob 5 % 95 % 25 % 75 % Pr(>|z|) stderr 
Andulay     0.333 0.091 0.663 0.214 0.469 0.42349216 0.192 
Antulang     0.333 0.112 0.888 0.304 0.696 1.00000000 0.192 
Basak      0.833 0.548 0.993 0.802 0.964 0.09916496 0.152 
Dauin Poblacion District 1 1.000 0.000 NA 0.000 1.000 0.99097988 0.000 
Guinsuan     0.500 0.223 0.940 0.474 0.819 0.56032414 0.204 
Kookoo's Nest    0.833 0.548 0.993 0.802 0.964 0.09916496 0.152 
Lutoban Pier    0.167 0.028 0.788 0.130 0.501 0.51171512 0.152 
Lutoban South    0.167 0.028 0.788 0.130 0.501 0.51171512 0.152 
Malatapay Pier    0.667 0.364 0.972 0.640 0.903 0.25767454 0.192

所以我的問題是雙重的：

是概率和置信區間正確的計算？
如何在bloxplot（盒子和鬍鬚圖）中繪製此圖？

EDIT這裏是經由dput一些示例數據（其還修改上面的表，以匹配數據）：

# dput(data[c("y", "site_name")]) 
data <- structure(list(y = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L), .Label = c("absent", "present"), class = "factor"), site_name = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 9L, 9L, 9L, 9L, 9L, 9L, 4L, 4L, 4L, 4L, 4L, 4L, 8L, 8L, 8L, 8L, 8L, 8L, 7L, 7L, 7L, 7L, 7L, 7L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("Andulay", "Antulang", "Basak", "Dauin Poblacion District 1", "Guinsuan", "Kookoo's Nest", "Lutoban Pier", "Lutoban South", "Malatapay Pier"), class = "factor")), .Names = c("y", "site_name"), row.names = c(125L, 123L, 126L, 124L, 128L, 127L, 154L, 159L, 157L, 158L, 156L, 155L, 111L, 114L, 116L, 115L, 112L, 113L, 152L, 151L, 148L, 150L, 153L, 149L, 143L, 146L, 144L, 147L, 142L, 145L, 164L, 165L, 161L, 163L, 160L, 162L, 120L, 122L, 121L, 117L, 118L, 119L, 137L, 136L, 139L, 141L, 140L, 138L, 129L, 134L, 131L, 135L, 133L, 130L), class = "data.frame") 
#

來源

2016-02-06 Dolf Andringa

任何可重複的例子......？ –

我添加了一些數據。我不知道我在上面的表格中列出了哪些種類，所以我更新了它們以匹配粘貼的數據。但問題本質上更具技術性。我如何獲得邏輯迴歸的預測值和置信區間的箱線圖？除非我理解完全錯誤的東西（不可想象） –

嗨，爲了提供您的數據，您可以用'dput（data [c（「y」，「site_name」）]）的結果來編輯您的問題。（希望人們可以將您的數據從您的問題複製到他們的R會話 - 我們無法使用您發佈的格式進行此操作） – user20650

這是一個僅基包最低公分母，，解決方案。

擬合模型：

mm <- glm(y~site_name,data=dd,family=binomial)

化妝與網站名稱的預測幀：

pframe <- data.frame(site_name=unique(dd$site_name))

預測（在分對數/線性預測器規模）中，用標準誤差

pp <- predict(mm,newdata=pframe,se.fit=TRUE) 
linkinv <- family(mm)$linkinv ## inverse-link function

將預測，下限和上限以及迴歸轉換爲概率尺度：

pframe$pred0 <- pp$fit 
pframe$pred <- linkinv(pp$fit) 
alpha <- 0.95 
sc <- abs(qnorm((1-alpha)/2)) ## Normal approx. to likelihood 
alpha2 <- 0.5 
sc2 <- abs(qnorm((1-alpha2)/2)) ## Normal approx. to likelihood 
pframe <- transform(pframe, 
        lwr=linkinv(pred0-sc*pp$se.fit), 
        upr=linkinv(pred0+sc*pp$se.fit), 
        lwr2=linkinv(pred0-sc2*pp$se.fit), 
        upr2=linkinv(pred0+sc2*pp$se.fit))

劇情。

with(pframe, 
{ 
    plot(site_name,pred,ylim=c(0,1)) 
    arrows(as.numeric(site_name),lwr,as.numeric(site_name),upr, 
      angle=90,code=3,length=0.1) 
})

由於箱線圖：

with(pframe, 
{ 
    bxp(list(stats=rbind(lwr,lwr2,pred,upr2,upr), 
      n = rep(1,nrow(pframe)), 
      conf = NA, 
      out = NULL, 
      group = NULL, 
      names=as.character(site_name))) 
})

有許多其他方法可以做到這一點;我建議

library("ggplot2") 
ggplot(pframe,aes(site_name,pred))+ 
    geom_pointrange(aes(ymin=lwr,ymax=upr))+ 
    geom_linerange(aes(ymin=lwr2,ymax=upr2),lwd=1.5)+ 
    coord_flip()

另一種解決方案是通過y~site_name-1，在這種情況下將要分配給每個站點的概率的單獨的參數，以適應模型，並使用profile()/confint()找到置信區間;這比依靠上述答案中參數/預測的抽樣分佈的正態性稍微更準確一些。

來源

2016-02-08 04:48:09

如何繪製物流GLM預測的R值和置信區間

回答

相關問題