(sigmoid) curve fitting glm in r

Question

I wish to visualize the relationship between my response variable, detection probability (P.det) and predictor variable (distance) for two categories (transmitter), show error bars and draw a (sigmoidal) curve through the averaged data points.

The dataset is like this:

df <- structure(list(distance = c(50L, 100L, 200L, 300L, 400L, 50L, 
100L, 200L, 300L, 400L), Transmitter = structure(c(1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L), .Label = c("CT", "PT"), class = "factor"), 
    P.det = c(0.918209097, 0.88375438, 0.709288774, 0.534977488, 
    0.341724516, 0.828123952, 0.822201191, 0.543289433, 0.352886247, 
    0.10082457), st.error = c(0.01261614, 0.014990469, 0.024136478, 
    0.027311169, 0.026941438, 0.018663591, 0.019420587, 0.02754911, 
    0.026809247, 0.017041264), ly = c(0.905592958, 0.868763911, 
    0.685152295, 0.50766632, 0.314783078, 0.809460361, 0.802780604, 
    0.515740323, 0.326077, 0.083783306), uy = c(0.930825237, 
    0.898744849, 0.733425252, 0.562288657, 0.368665955, 0.846787544, 
    0.841621778, 0.570838544, 0.379695494, 0.117865833), Valid.detections = c(18, 
    12.5472973, 8.608108108, 4.287162162, 2.158783784, 12.46959459, 
    7.956081081, 4.550675676, 1.682432432, 0.39527027), False.detections = c(0.388513514, 
    0.550675676, 0.368243243, 0.263513514, 0.131756757, 0.533783784, 
    0.385135135, 0.277027027, 0.182432432, 0.14527027)), .Names = c("distance", 
"Transmitter", "P.det", "st.error", "ly", "uy", "Valid.detections", 
"False.detections"), class = "data.frame", row.names = c(NA, 
-10L))

I managed to get the first 2 parts done, but am stuck at the last part. The code to draw the graph with error bars:

library(lattice)
    prepanel.ci <- function(x, y, ly, uy, subscripts, ...)
    {
      x <- as.numeric(x)
      ly <- as.numeric(ly[subscripts])
      uy <- as.numeric(uy[subscripts])
      list(ylim = range(y, uy, ly, finite = TRUE))
    }

    panel.ci <- function(x, y, ly, uy, subscripts, pch = 16, ...)
    {
      x <- as.numeric(x)
      y <- as.numeric(y)
      ly <- as.numeric(ly[subscripts])
      uy <- as.numeric(uy[subscripts])
      panel.arrows(x, ly, x, uy, col = "black",
                   length = 0.25, unit = "native",
                   angle = 90, code = 3)
      panel.xyplot(x, y, pch = pch, ...)
    }

xyplot(P.det~distance, type=c("p","g"),
       ylim=c(0,1),
       ylab="Detection probability", xlab="Distance (m)", 
       group=Transmitter,
       data=df,
       ly = df$ly,
       uy = df$uy,
       prepanel = prepanel.ci,
       panel = panel.superpose,
       panel.groups = panel.ci,
       col=c(1,1),
       layout=c(1,1),
       between=list(x=2),
       scales=list(x=list(alternating=c(1,1), tck=c(1,0)),y=list(alternating=c(1,1), tck=c(1,0))), # ticks inside = tck=c(-1,0)
       aspect=1,
       main="Detection probability vs distance per transmitter type",
)

The reason why I state "glm" in the title is because the data analysis was carried out using a binomial glm() using the lme4 package.

I noticed another thread which is similar to mine: find the intersection of abline with fitted curve, however the difference is that while my graph is also based on 1 y per 1 x, my glm is based on a multitude of y's per x. So following the same codes in this thread returns an error stating that the lengths are not of equal length. It also doesn't seem to work for an "xyplot".

Thanks

score 4 · Accepted Answer · edited May 23 '17 at 12:07

4

This is fairly straightforward using ggplot:

library(ggplot2)
ggplot(data = df, aes(x = distance, y = P.det, colour = Transmitter)) +
  geom_pointrange(aes(ymin = P.det - st.error, ymax = P.det + st.error)) +
  geom_smooth(method = "glm", family = binomial, se = FALSE)

enter image description here

Regarding the glmwarning message, see e.g. here.

edited May 23 '17 at 12:07

Community

1
1

answered Jul 09 '14 at 14:59

Henrik

65,555
14
143
159

Thanks, do you also know how to do this with xyplot? I would like to use it in my publication, and actually ggplot doesn't seem to lend itself for this. – FlyingDutch Jul 09 '14 at 15:20
Sorry, I am not that familiar with `lattice`. But it is not difficult to customize your `ggplot`. You may have a look at [**`theme`**](http://docs.ggplot2.org/current/theme.html) and e.g. [**`theme_bw`**](http://docs.ggplot2.org/current/theme_bw.html), `scale_colour_xxx`, `scale_linetype_xxx`, `scale_shape_xxx` etc. See also [**"How to format plots for publication using ggplot2"**](http://www.noamross.net/blog/2013/11/20/formatting-plots-for-pubs.html) – Henrik Jul 09 '14 at 15:20
Right, I am only starting to know lattice, I guess I'm wrong about ggplot. I'll mark your answer as accepted and look into ggplot more. Thanks! – FlyingDutch Jul 09 '14 at 15:21
Other good tutorials for ggplot: [**here**](http://www.cookbook-r.com/Graphs/) and [**here**](http://sape.inf.usi.ch/quick-reference/ggplot2). – Henrik Jul 09 '14 at 15:26

score 1 · Answer 2 · answered Jul 09 '14 at 15:53

In Lattice, you might add the smoothing into your custom panel function. You can change it to include

panel.ci <- function(x, y, ly, uy, subscripts, type="p", pch = 16, ...)
{
  x <- as.numeric(x)
  y <- as.numeric(y)
  ly <- as.numeric(ly[subscripts])
  uy <- as.numeric(uy[subscripts])
  panel.arrows(x, ly, x, uy, col = "black",
               length = 0.25, unit = "native",
               angle = 90, code = 3)
  panel.xyplot(x, y, pch = pch, type=type, ...)

  #calculate smooth curve
  gg <- glm(y~x, family="binomial")
  gx <- seq(min(x), max(x), length.out=100)
  panel.xyplot(gx, predict(gg, data.frame(x=gx), type="response"),
      pch = pch, type="l", ...)
}

Here we do the glm ourselves and draw the response on the plot.

enter image description here

(sigmoid) curve fitting glm in r

2 Answers2

Linked