R Finding logistic curve with nls

Question

I have a problem with logistic curve in R for panel data which are here:

https://docs.google.com/spreadsheets/d/1SO3EzFib7T3XqTz1xZCU2bZFTB0Ddhou/edit?usp=sharing&ouid=110784858039906954607&rtpof=true&sd=true

I tried:

log <- nls(lrc~SSlogis(time, Asym, xmid, scal), data = data_log)

I have an error: 'qr.solve(QR.B, cc)':singular matrix 'a' in solve. What can I do?

Answer 1

I get a different error:

Error in nls(y ~ 1/(1 + exp((xmid - x)/scal)), data = xy, start = list(xmid = aux[[1L]], : step factor 0.000488281 reduced below 'minFactor' of 0.000976562

(it's not surprising that different platforms will get slightly different results on numerically "difficult" problems ...)

However, here's what plot(lrc ~ time, data = dd) produces when I use your data:

It seems optimistic to think that you could fit a logistic curve to these data, or that the fit would make very much sense ...

I did find that I could fit a logistic to the logged data, i.e. nls(log(lrc) ...)

plot(log(lrc) ~ time, data =dd)
tvec <- seq(2008, 2016, by = 0.1)
lines(tvec, predict(m2, newdata=data.frame(time=tvec)), col=2, lwd=2)

If I really needed logistic coefficients for this plot (e.g. to compare with other cases) I would fit a linear model to the data, assume that the midpoint was equal to the midpoint of the data, set the scaling parameter equal to the linear slope coefficient divided by 4 (this is a standard rule of thumb for the logistic), and say that the asymptotic value could not be estimated.

---

Plotting all of your data unit-by-unit doesn't make much more hopeful that you will be able to fit logistic curves unit-by-unit either. While there might be a few units where a logistic curve is a sensible description of the data, it's not in most cases. You might have to back up and consider your analytical strategy — i.e., what are you hoping to learn from these data, and how might you go about it? (If you can frame the question suitably you could post on CrossValidated. If you are a student/trainee, you might want to ask your supervisor/teacher/mentor for advice ...

library(readxl)
library(tidyverse)
(dd <- read_excel("data2.xlsx")
  %>% pivot_longer(-code, names_transform = as.numeric,
                   names_to = "year")
  ## indices to break groups into chunks/facets
  %>% mutate(grp_cat = factor(as.numeric(factor(code)) %% 30))
)
gg1 <- ggplot(dd, aes(year, value, group = code,
               colour=code)) + geom_point() +
  facet_wrap(~grp_cat, scale="free_y") +
  expand_limits(y=0) +
  theme_bw() +
  theme(legend.position = "none",
        panel.spacing = grid::unit(0, "lines"),
        axis.text.x = element_blank())
gg1 + geom_smooth(se=FALSE)

ggsave("all.png", width=8, height=8)