4
$\begingroup$

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

enter image description here

$\endgroup$

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.