0
$\begingroup$

I seem to occasionally find datasets where the relationship between X and Y is concave downward. It seems like it should be trivial to find a link function that fits a concave downward curve, but they all seem to fit concave upward curves. Here's a reproducible example (assuming there's no problems downloading my flexplot package):

#devtools::install_github('dustinfife/flexplot')
require(flexplot)
require(tidyverse)

        ### read in data
d <- read.csv("https://raw.githubusercontent.com/dustinfife/flexplot/master/data/nsduh_drug_users.csv")


        ### plot raw data with lowess line. Curve is concave downward
a <- flexplot::flexplot(distress~MI, data=d)

        ### fit a model
mod <- glm(distress~MI, data=d %>% mutate(distress = distress+1), family=Gamma(link="log"))
        ### visualize the fit of the model (Curve is concave upward)
b <- flexplot::compare.fits(distress~MI, data=d, model1=mod)
cowplot::plot_grid(a,b)

And here are the results of the two plots:

concave downward vs upward plot

Notice that the left loess line shows a concave downward curve, while the right plot shows a concave upward curve.

I looked into several options, including:

  1. A logistic growth curve with a set carrying capacity. This seems to be the type of curve that I want and I can model it in nlm, but I figured I wanted to keep it within glm. Why? Because nlm fits using least squares and (I assume) it assumes a normal distribution of residuals. Judging from the plot, the residuals disperse on a log scale, not a normal scale.

  2. Set my own custom link function in glm. I don't know why, but that idea makes me cringe. Seems complicated, but if it's what's necessary, I'll do it.

  3. Transform the data to be in the 0-1 interval and fit some sort of a logistics. I actually tried this, but R barked at me because the data were not binary. But I assume there's some way of overcoming that, although I'm not sure if it's a good idea because the residuals would not be distributed as a logistic.

  4. (Ideal solution). Find some link function that already exists in some package and use that with a gamma distribution.

So, my questions are:

  • Is there a link function in glm (or some external package) that allows concave downward fits?
  • Is my thinking correct that nlm won't work because the residuals exhibit heteroskedasticity?
  • If no to #1, is there an easy way to fit a concave downward function that models the residuals appropriately?

Edit: To show an example that is more clearly concave downward:

This example is borrowed from this link, but I modified the dataset to add heteroscedasticity:

#Here's the data
mass<-c(6.25,10,20,23,26,27.6,29.8,31.6,37.2,41.2,48.7,54,54,63,66,72,72.2,76,75) #Wilson's mass in pounds
days.since.birth<-c(31,62,93,99,107,113,121,127,148,161,180,214,221,307,452,482,923, 955,1308) #days since Wilson's birth
data<-data.frame(mass,days.since.birth) #create the data frame

    #### add heteroscedasticity, conditional on mass
mass <- mass + rnorm(length(mass), 0, 10*(mass/max(mass)))  

    #### fit with nlme
wilson<-nls(mass~phi1/(1+exp(-(phi2+phi3*days.since.birth))), 
        start=list(phi1=100,phi2=-1.096,phi3=.002),data=data,trace=TRUE)

    ### prepare for ggplot
phi1<-coef(wilson)[1]
phi2<-coef(wilson)[2]
phi3<-coef(wilson)[3]
x<-c(min(data$days.since.birth):max(data$days.since.birth)) #construct a range of x values bounded by the data
y<-phi1/(1+exp(-(phi2+phi3*x))) #predicted mass
predict<-data.frame(x,y) #create the prediction data frame#And add a nice plot (I cheated and added the awesome inset jpg in another program)
a <- ggplot(data=data,aes(x=days.since.birth,y=mass))+
    geom_point()+
    labs(x='Days Since Birth',y='Mass (lbs)')+
    scale_x_continuous(breaks=c(0,250,500,750, 1000,1250))+
    scale_y_continuous(breaks=c(0,10,20,30,40,50,60,70,80))+
    geom_line(data=predict,aes(x=x,y=y), size=1) +
    theme_bw()

        #### fit with Gamma
mod <- glm(mass~days.since.birth, data=data, family=Gamma(link="log"))
b <- compare.fits(mass~days.since.birth, data=data, model1=mod)
cowplot::plot_grid(a,b)

And the results:

logistic growth model vs. gamma

This is a gamma family with a log link, but I have also tried other link functions (1/mu^2, inverse, logit, log) and they all are either concave upward or R throws an error message (because it's not within a 0-1 interval).

$\endgroup$
  • $\begingroup$ Given the data you present, it's a little difficult to be convinced that the fitted curve should be either concave or convex. Can you replicate this problem with more-clear data? $\endgroup$ – Sal Mangiafico Sep 15 '19 at 13:38
  • $\begingroup$ With the additional example, it's an interesting problem.... If I have time, I'll play with it a bit.... $\endgroup$ – Sal Mangiafico Sep 16 '19 at 0:03
  • $\begingroup$ It's unclear what you mean by "concave downward" given that you characterize curves that do have concave portions, such as the logit and log, as not being concave downward, and you include both the log and (apparently) the reciprocal square as not being concave downward when these have opposite concavities. Please tell us what you mean by this term. $\endgroup$ – whuber Sep 16 '19 at 10:59
  • $\begingroup$ @whuber - I just want the curve to rapidly increase then taper off as it approaches some asymptote, like the left image in the second figure. Any glm curve I try to fit does the opposite: it increases slowly as X increases, then accelerates it's growth. Does that help? $\endgroup$ – dfife Sep 16 '19 at 12:58
  • $\begingroup$ It doesn't help, because the fact remains that glm already supports plenty of links like you describe. Apparently, then, something is going wrong. A glance at your last code block shows you aren't even invoking a GLM: you're using nls, which is quite a different procedure! The problem probably isn't the link at all, but lies rather in the inflexible model you are using: spline the data. $\endgroup$ – whuber Sep 16 '19 at 13:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.