Cannot plot p-value on simple logistic regression

Question

I am trying to plot a simple logistic regression in R.

I am following this tutorial to conduct the logistic regression and calculate a P-value (https://mgimond.github.io/Stats-in-R/Logistic.html). I am trying to use ggplot2 and ggpmisc to plot the regression. I have been trying to use this guide (http://cran.nexr.com/web/packages/ggpmisc/vignettes/user-guide-1.html#stat_fit_glance) to stat_fit_glance to display a p-value

require(cowplot)
require(ggplot2)
library(ggpmisc)
library(rms)

dataset=read.table('input.txt', header=TRUE)
model <- glm(variable ~ ancestry, data=dataset, family=binomial)
summary(model)
#plot logistic regression curve
plot <- ggplot(dataset, aes(x=ancestry, y=variable)) + 
  geom_point(alpha=.5, color=dataset$colorsite) +
  stat_smooth(method="glm", se=FALSE, method.args = list(family=binomial)) + stat_fit_glance(method = "glm", method.args = list(formula = formula), geom = "text", aes(label = paste("P-value = ", signif(..p.value.., digits = 4), sep = ""))) 
ggsave("output.pdf")

The output however comes out as

> source("C:/Users/Deven/Desktop/logistic/script.R")
Saving 7 x 7 in image
`geom_smooth()` using formula 'y ~ x'
Warning message:
Computation failed in `stat_fit_glance()`:
object of type 'closure' is not subsettable

I have also tried stat_cor from ggpubr, but that seem to be generating different p-values and R^2 values from what I have calculated.

UPDATE BASED ON COMMENTS: + stat_poly_eq(formula = y ~ x, method="glm", aes(x = ancestry, y = variable, label = paste(..p.value.label..,sep = "~~~~")),parse = TRUE) fails due to

1: Computation failed in `stat_poly_eq()`:
Method 'glm' not yet implemented. 

If I remove method it defaults to a linear regresssion (and gives p values that do not correspond to a logistic regression).

SECOND UPDATE

model <- glm(variable ~ ancestry, data=dataset, family=binomial)
summary(model)
#plot logistic regression curve
plot <- ggplot(dataset, aes(x=ancestry, y=variable)) +
  geom_point(alpha=.5, color=dataset$colorsite) +
  stat_smooth(method="glm", se=FALSE, method.args = list(family=binomial)) + stat_fit_tidy(method = "glm",method.args = list(family=binomial,formula=y~x), mapping = aes(label = sprintf("Coef = %.3g\np-value = %.3g",after_stat(x_estimate),after_stat(x_p.value))))
ggsave("variable.pdf")

yields the following error:

Saving 7 x 7 in image
`geom_smooth()` using formula 'y ~ x'
Warning message:
Computation failed in `stat_fit_tidy()`:
no applicable method for 'tidy' applied to an object of class "c('glm', 'lm')" 

YET ANOTHER UPDATE

library(ggplot2)
library(ggpmisc)

da =read.table('data.txt', header=TRUE)


model = glm(variable ~ ancestry,family=binomial,data=da)
summary(model)

ggplot(da,aes(x = ancestry,y = variable)) + geom_point() +
  stat_smooth(method="glm",se=FALSE,method.args = list(family=binomial)) +
  stat_fit_tidy(method = "glm",method.args = list(family=binomial,formula=y~x),
                mapping = aes(label = sprintf("Coef = %.3g\np-value = %.3g",
                                              after_stat(x_estimate),after_stat(x_p.value))))
ggsave("test.pdf")

works in theory, but the p-value it gives me is very different from the p value that I calculated manually (which corresponds to the one from lrm(variable ~ ancestry, dataset))...

Not sure at all what is going on here...

Answer 1

Thank you for all the help, but unfortunately nothing automated worked, so I came up with this instead

require(cowplot)
require(ggplot2)
library(ggpmisc)
library(rms)


dataset=read.table('data.txt', header=TRUE)

model <- glm(variable ~ ancestry, data=dataset, family=binomial)
summary(model)
M1 <- glm(variable ~ ancestry, dataset, family = binomial)
M1
M1$null.deviance
M1$deviance
modelChi <- M1$null.deviance - M1$deviance
pseudo.R2 <- modelChi / M1$null.deviance
pseudo.R2
test <-lrm(variable ~ ancestry, dataset)
Chidf <- M1$df.null - M1$df.residual
chisq.prob <- 1 - pchisq(modelChi, Chidf)
chisq.prob

#plot logistic regression curve
all_variable <- ggplot(dataset, aes(x=ancestry, y=variable)) +
  geom_point(alpha=.5, color=dataset$colorsite) +
  stat_smooth(method="glm", se=FALSE, method.args = list(family=binomial)) + annotate("text", x=-Inf, y=Inf, hjust = 0, vjust = 2.5, label=paste("p-value: ",signif(chisq.prob, digits = 3),"\nR2: ",signif(pseudo.R2, digits = 3),sep="") )+
  ggtitle("Title not relevant to Stack Overflow") 
ggsave("variable.pdf")

Answer 2

There is a table on ggpmisc help page that specifies what can be applied to each type of models.

You have a glm, so glance() from tidy will not give you a p-value. Using an example:

library(ggplot2)
library(ggpmisc)

da = MASS::Pima.tr
da$label = as.numeric(da$type=="Yes")

model = glm(label ~ bmi,family=binomial,data=da)
summary(model)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -4.11156    0.92806  -4.430 9.41e-06 ***
bmi          0.10482    0.02738   3.829 0.000129 ***

You can see glance will not give you a p-value :

broom::glance(model)
# A tibble: 1 x 8
  null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1          256.     199  -120.  244.  251.     240.         198   200

You need to use tidy() and as @JonSpring mentioned in the comment, provide the formula, so something like this:

ggplot(da,aes(x = bmi,y = label)) + geom_point() +
stat_smooth(method="glm",se=FALSE,method.args = list(family=binomial)) +
stat_fit_tidy(method = "glm",method.args = list(family=binomial,formula=y~x),
mapping = aes(label = sprintf("Coef = %.3g\np-value = %.3g",
after_stat(x_estimate),after_stat(x_p.value))))