2017-03-12 25 views
0

使用R,我想輸入如何找到預測邏輯曲線的方程(從經驗數據預測)。如何找到預測對數曲線的方程?

我能找到的最接近的命令是library(ggpmisc)stat_poly_eq函數,並在使用ggplot2繪製曲線時使用。但是,我只能用這種方法在圖上打印y = 0.48,但我需要整個方程。

下面是一些供參考的R代碼裏面:

gg.disp.adults<- ggplot(sub.data, aes(x=SVL3, y=Disp01)) + 
    geom_point(size=4) + 
    stat_smooth(aes(y= Disp01), method="glm", method.args=list(family="binomial"), se=F) + 
    stat_poly_eq(aes(label=paste(..eq.label..,..rr.label..,sep="~~~~")), 
      rr.digits=3, coef.digits=2, 
      formula = y~1/(1+exp(-x)), 
      parse = TRUE) + 
theme(axis.text.x=element_text(size=14, color="black"), 
    axis.text.y=element_text(size=14, color="black"), 
    axis.line=element_line(size=1), 
    axis.title.x=element_text(size=14), 
    axis.title.y=element_text(size=14), 
    panel.background=element_rect(fill="white")) + 
ylab("Dispersal Probability") + 
ylim(0,1)+ 
xlab("Adult SVL") 
gg.disp.adults 

plot

請如何去獲得實際的方程預測曲線建議。

編輯:

我能夠找到使用係數從輸出摘要的預測曲線的方程,並將其應用到一個標準的邏輯斯諦方程: Y〜1 /(1 + EXP(-x ))。以下是R代碼供參考:

Call: 
glm(formula = Disp01 ~ SVL3, family = binomial, data = sub.data) 

Deviance Residuals: 
    Min  1Q Median  3Q  Max 
-2.1207 -0.8812 -0.4844 0.8885 1.9168 

Coefficients: 
      Estimate Std. Error z value Pr(>|z|)  
(Intercept) -12.1800  3.6968 -3.295 0.000985 *** 
SVL3   0.1845  0.0561 3.289 0.001006 ** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

x<- seq(50,80,1) 
y<- 1/(1+exp(-(0.1845*x-12.18))) #coefficients from glm2 output 
plot(y~x, typ="l", ylim=c(0,1)) 
points(Disp01~SVL3, data=sub.data, pch=16, cex=1.5) 

...我無法發佈新預測曲線的照片。

回答

1

由於stat_smooth只是對數據運行glm(y ~ x, ...),你可以做恢復模型估計值相同。下面是一些由數據爲例:

化妝數據:

library(ggplot2) 
set.seed(123) 
d <- data.frame(x = rnorm(100, 0, 2)) 
d$y <- rbinom(100, 1, plogis(d$x)) 

適合使用glm。因此fit會報告你想要的公式:

fit <- glm(y ~ x, data = d, family = binomial(link = "logit")) 
fitted <- data.frame(x = d$x, y = predict(fit, type = "response")) 

在這裏我們可以看到,預測值從fitted搭配起來正好與stat_smooth行:

ggplot(d, aes(x, y)) + 
    stat_smooth(method = "glm", method.args = list(family = "binomial")) + 
    geom_point(data = fitted, aes(x, y)) 

enter image description here