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")
