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
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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
# pint.R: plotmo functions for confidence and prediction intervals
# Handle plotmo's "level" argument. Return a prediction interval dataframe
# with either or both of the following sets of columns. What columns get
# returned depends on the capabilities of the object's predict method.
# For example, predict.lm allows us to return both i and ii, and for
# earth models we can return only i.
#
# (i) lwr, upr intervals for prediction of new data
#
# (ii) cint.lwr, cint.upr intervals for prediction of mean response
plotmo_pint <- function(object, newdata, type, level, trace, ipred, inverse.func)
{
if(!is.specified(level))
return(NULL)
trace2(trace, "plotmo_pint for %s object\n", class.as.char(object))
stopifnot.string(type)
# call plotmo.pint.xxx where xxx is object's class
intervals <- plotmo.pint(object, newdata, type, level, trace)
if(!is.null(intervals$lwr)) {
intervals$lwr <-
apply.inverse.func(inverse.func, intervals$lwr, object, trace)
intervals$upr <-
apply.inverse.func(inverse.func, intervals$upr, object, trace)
}
if(!is.null(intervals$cint.lwr)) {
intervals$cint.lwr <-
apply.inverse.func(inverse.func, intervals$cint.lwr, object, trace)
intervals$cint.upr <-
apply.inverse.func(inverse.func, intervals$cint.upr, object, trace)
}
print_summary(intervals, "prediction intervals", trace)
intervals
}
# Return a data.frame with either or both of the following variables:
# (i) lwr, upr intervals for prediction of new data
# (ii) cint.lwr, cint.upr intervals for prediction of mean response
plotmo.pint <- function(object, newdata, type, level, trace, ...)
{
UseMethod("plotmo.pint")
}
plotmo.pint.default <- function(object, newdata, type, level, trace, ...)
{
stop0("the level argument is not supported for ",
class.as.char(object, quotify=TRUE), " objects")
}
plotmo.pint.lm <- function(object, newdata, type, level, trace, ...)
{
# lm objects with weights do not support confidence intervals on new data
if(!is.null(object$weights))
stop0("the level argument is not supported on lm objects with weights")
pints <- predict(object, newdata, interval="prediction", level=level)
cints <- predict(object, newdata, interval="confidence", level=level)
data.frame(
lwr = pints[,"lwr"], # intervals for prediction of new data
upr = pints[,"upr"],
cint.lwr = cints[,"lwr"], # intervals for prediction of mean response
cint.upr = cints[,"upr"])
}
plotmo.pint.glm <- function(object, newdata, type, level, trace, ...)
{
predict <- try(predict(object, newdata, type="link", se.fit=TRUE))
if(is.try.err(predict))
stopf("\"predict(mod, newdata, type=\"link\", se.fit=TRUE)\" failed for \"%s\" model", class(object)[1])
ilink = try(family(object)$linkinv)
if(is.try.err(ilink))
stopf("\"family(mod)$linkinv\" failed for \"%s\" model", class(object)[1])
df <- try(df.residual(object))
if(is.try.err(df))
stopf("\"df.residual(mod)\" failed for \"%s\" model", class(object)[1])
q <- get.quant(level, df=df, trace=trace>=2) # e.g for level=.95 return 1.96 if df=Inf
if(!is.null(object$weights) && !all(object$weights == object$weights[1]))
warnf("the level argument may not work correctly on glm objects built with weights")
data.frame(cint.lwr = ilink(predict$fit - q * predict$se.fit), # changed in plotmo 3.7.0, jan 2026
cint.upr = ilink(predict$fit + q * predict$se.fit))
}
# package mgcv, or package gam version less than 1.15
plotmo.pint.gam <- function(object, newdata, type, level, trace, ...)
{
quant <- 1 - (1 - level) / 2 # .95 becomes .975
predict <- predict(object, newdata, type=type, se.fit=TRUE)
# special handling for where user used gam::gam instead of mgcv::gam
# (applies only package gam version less than 1.15)
if(class(predict)[1] == "numeric" && "package:gam" %in% search()) {
cat("\n")
stop0("gam objects in the 'gam' package do not support ",
"confidence intervals on new data")
}
if(!is.null(object$weights) && !all(object$weights == object$weights[1]))
warnf(
"the level argument may not work correctly on gam objects built with weights")
data.frame(cint.lwr = predict$fit - quant * predict$se.fit,
cint.upr = predict$fit + quant * predict$se.fit)
}
# package gam version 1.15 or higher
plotmo.pint.Gam <- function(object, newdata, type, level, trace, ...)
{
quant <- 1 - (1 - level) / 2 # .95 becomes .975
predict <- predict(object, newdata, type=type, se.fit=TRUE)
if(class(predict)[1] == "numeric") {
cat("\n")
stop0("Gam objects do not support confidence intervals on new data")
}
if(!is.null(object$weights) && !all(object$weights == object$weights[1]))
warnf(
"the level argument may not work correctly on Gam objects built with weights")
data.frame(cint.lwr = predict$fit - quant * predict$se.fit,
cint.upr = predict$fit + quant * predict$se.fit)
}
plotmo.pint.quantregForest <- function(object, newdata, type, level, trace, ...)
{
q0 <- (1 - level) / 2 # .95 becomes .025
q1 <- 1 - q0 # .975
predict <- predict(object, newdata, quantiles=c(q0, q1))
data.frame(lwr = predict[,1], upr = predict[,2])
}
plotmo.pint.earth <- function(object, newdata, type, level, trace, ...)
{
pints <- predict(object, newdata=newdata, type=type, interval="pint", level=level)
if(is.null(newdata)) {
cints <- predict(object, newdata=NULL, type=type, interval="cint", level=level)
pints$cint.upr <- cints$upr
pints$cint.lwr <- cints$lwr
}
pints
}
|
