# How to calculate standard errors of a non-linear model prediction?

I'm trying to understand how to show the prediction error of a model fit in R using the non-linear least squares function nls. Although there is an argument se.fit=TRUE in thepredict.nls function, the help explains "At present se.fit and interval are ignored." I'm able to calculate the confidence interval of the fitted parameters, but am unsure if this information can be used to calculate prediction std. error. Below is an example:

### Example data:

df <- structure(list(L = c(13L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,
15L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
16L, 16L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L,
17L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L,
19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,
20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L,
21L, 21L, 21L, 21L, 21L, 22L, 22L, 22L, 22L, 22L, 22L, 23L, 23L,
23L, 23L, 23L, 23L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L, 24L,
24L, 24L, 24L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 26L,
26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 26L, 27L,
27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L, 27L,
28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L,
28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L, 28L,
28L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L,
29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L, 29L,
30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L, 30L,
30L, 30L, 30L, 30L, 30L, 30L, 31L, 31L, 31L, 31L, 31L, 31L, 31L,
31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 31L, 32L,
32L, 32L, 32L, 32L, 32L, 33L, 33L, 34L, 34L, 34L, 34L, 34L, 34L,
34L, 35L, 35L, 35L, 35L, 37L, 38L), mat = c(0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L,
0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L,
1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L)), .Names = c("L", "mat"), class = "data.frame", row.names = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
"14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24",
"25", "26", "114", "115", "27", "28", "29", "30", "31", "32",
"33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43",
"44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54",
"55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65",
"116", "117", "66", "67", "68", "69", "70", "71", "72", "118",
"119", "73", "74", "75", "76", "120", "121", "77", "78", "79",
"122", "123", "124", "80", "81", "82", "83", "84", "85", "86",
"125", "126", "127", "128", "129", "87", "88", "89", "90", "130",
"131", "132", "133", "134", "91", "92", "93", "94", "95", "96",
"135", "136", "137", "138", "139", "140", "141", "97", "98",
"99", "100", "142", "143", "144", "145", "146", "147", "148",
"149", "150", "151", "101", "102", "103", "104", "105", "106",
"152", "153", "154", "155", "156", "157", "158", "159", "160",
"161", "162", "163", "164", "165", "166", "167", "168", "169",
"170", "171", "172", "107", "108", "109", "110", "111", "173",
"174", "175", "176", "177", "178", "179", "180", "181", "182",
"183", "184", "185", "186", "187", "188", "189", "190", "191",
"192", "112", "193", "194", "195", "196", "197", "198", "199",
"200", "201", "202", "203", "204", "205", "206", "207", "208",
"209", "210", "113", "211", "212", "213", "214", "215", "216",
"217", "218", "219", "220", "221", "222", "223", "224", "225",
"226", "227", "228", "229", "230", "231", "232", "233", "234",
"235", "236", "237", "238", "239", "240", "241", "242", "243",
"244", "245", "246", "247", "248", "249"))

df2 <- aggregate(df, by=list(as.factor(df$L)), FUN="mean")[,-1] names(df2) <- c("L", "pmat")  ### Fit the NLS model fmla <- "mat ~ 1 / (1 + exp(-a*(L-L50)))" fit <- nls( fmla, data=df, start=c(a=0.3, L50=25) ) summary(fit) fit.ci <- confint(fit) fit.ci newdat <- data.frame(L=seq(0,50,.1)) pred <- predict(fit, newdat)  Below, I compare this fit to a GLM model using a binomial distribution. predict.glm allows me to easily calculate the std. error of the prediction. You can see that the predictions between NLS and GLM are quite similar on this data. ### Comparison to binomial GLM fit2 <- glm(mat ~ L, data = df, family=binomial) pred2 <- predict(fit2, newdat, type="response", se.fit=TRUE) plot(df, xlim=c(0,50)) points(pmat ~ L, df2, pch=20, col=8) lines(newdat$L, pred)
lines(newdat$L, pred2$fit, col=2)
polygon(x=c(newdat$L, rev(newdat$L)), y=c(pred2$fit+pred2$se.fit, rev(pred2$fit-pred2$se.fit)), col=rgb(1,0,0,0.2), border=NA)
legend("right", legend=c("NLS", "binom. GLM"), col=1:2, lty=1, bty="n")
legend("left", legend=c("raw", "binned"), col=c(1,8), pch=c(1,20), bty="n")
`