View source: R/predict.gravity.R
predict.gravity | R Documentation |
predict method for class "gravity"
## S3 method for class 'gravity'
predict(
object,
newdata,
groups = NULL,
back.transform = c("none", "simple", "Miller", "Naihua"),
...
)
object |
Object of class gravity |
newdata |
New data used for obtaining the predictions, can be a data.frame or nffGroupedData |
groups |
Grouping factor acting as random effect. If used, must match levels used in model, otherwise leave it null and do not convert to groupedData |
back.transform |
Method to back transform data, default is none and log predictions will be returned. |
... |
Arguments passed to predict.lme or predict.lm |
Please note that the entire gravity equation is log transformed so, your parameter space is on a log scale, not just y. This means that for a meaningful prediction the "newdata" also needs to be on a log scale.
For the back.transform argument, the simple back-transform method uses the form exp(y-hat)0.5*variance whereas Miller uses exp(sigma)*0.5 as the multiplicative bias factor. Naihua regresses y~exp(y-hat) with no intercept and uses the resulting coefficient as the multiplicative bias factor. The Naihua method is intended for results with non-normal errors. You can check the functional form by simply plotting y (non-transformed) against the fit. The default is to output the log scaled predictions.
Vector of model predictions
Jeffrey S. Evans <jeffrey_evans@tnc.org> and Melanie A. Murphy <melanie.murphy@uwyo.edu>
Miller, D.M. (1984) Reducing Transformation Bias in Curve Fitting The American Statistician. 38(2):124-126
Naihua, D. (1983) Smearing Estimate: A Nonparametric Retransformation Method Journal of the American Statistical Association, 78(383):605–610.
library(nlme)
data(ralu.model)
back.transform <- function(y) exp(y + 0.5 * stats::var(y, na.rm=TRUE))
rmse = function(p, o){ sqrt(mean((p - o)^2)) }
x = c("DEPTH_F", "HLI_F", "CTI_F", "cti", "ffp")
sidx <- sample(1:nrow(ralu.model), 100)
train <- ralu.model[sidx,]
test <- ralu.model[-sidx,]
# Specify constrained gravity model
( gm <- gravity(y = "DPS", x = x, d = "DISTANCE", group = "FROM_SITE",
data = train, ln = FALSE) )
( p <- predict(gm, test[,c(x, "DISTANCE")]) )
rmse(back.transform(p), back.transform(ralu.model[,"DPS"][-sidx]))
# WIth model sigma-based back transformation
( p <- predict(gm, test[,c(x, "DISTANCE")], back.transform = "simple") )
( p <- predict(gm, test[,c(x, "DISTANCE")], back.transform = "Miller") )
( p <- predict(gm, test[,c(x, "DISTANCE")], back.transform = "Naihua") )
# Using grouped data
test <- nlme::groupedData(stats::as.formula(paste(paste("DPS", 1, sep = " ~ "),
"FROM_SITE", sep = " | ")),
data = test[,c("DPS", "FROM_SITE", x, "DISTANCE")])
( p <- predict(gm, test, groups = "FROM_SITE") )
( y.hat <- back.transform(ralu.model[,"DPS"][-sidx]) )
na.idx <- which(is.na(p))
rmse(back.transform(p)[-na.idx], y.hat[-na.idx])
# Specify unconstrained gravity model (generally, not recommended)
( gm <- gravity(y = "DPS", x = x, d = "DISTANCE", group = "FROM_SITE",
data = train, ln = FALSE, constrained=TRUE) )
( p <- predict(gm, test[,c(x, "DISTANCE")]) )
rmse(back.transform(p), back.transform(ralu.model[,"DPS"][-sidx]))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.