Search code examples
rregressionnlsextrapolation

Extrapolation in R

I have a reference data set, p, 

p=structure(list(v = 0:26, t = c(Inf, 1.016, 0.568, 0.418666666666667, 
0.344, 0.2992, 0.269333333333333, 0.248, 0.232, 0.219555555555556, 
0.2096, 0.201454545454545, 0.194666666666667, 0.188923076923077, 
0.184, 0.179733333333333, 0.176, 0.172705882352941, 0.169777777777778, 
0.167157894736842, 0.1648, 0.162666666666667, 0.160727272727273, 
0.15895652173913, 0.157333333333333, 0.15584, 0.154461538461538
)), .Names = c("v", "t"), row.names = c(NA, -27L), class = "data.frame")

Which p$v and p$t have the following relation :

t=(0.16*(0.75*v + 5.6))/v

my second data set is measured data ,w, contains the same variable like :

w=structure(list(v = c(0, 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), t = c(0.235291176470588, 
0.354020375722543, 0.310974343434343, 0.25272725498699, 0.20351968240702, 
0.163155804025208, 0.132330740162655, 0.108593108108108, 0.0859813015873016, 
0.0655131899302683, 0.0492580103144236, 0.0368029846567365, 0.030538003253355, 
0.0300744415648525, 0.0347586421891237, 0.0451097744360902, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)), .Names = c("v", "t"
), row.names = c(NA, -27L), class = "data.frame")

My reference data p, follows a power law, I would like to do a prediction based on power law fit on the data to replace the NA on my measured data. How could I do this in R ?

Solution

  • How about something like this?

    # Non-linear power law fit
fit.nls <- nls(
    t ~ a * v^b,
    data = w[-1, ],
    start = list(a = 0.4, b = -0.7),
    na.action = na.exclude);
# Linear fit with log-log transformation
fit.lm <- lm(
    log(t) ~ log(v),
    data = w[-1, ])

# Plot
w %>%
    rename(t.src = t) %>%
    mutate(
        t.pred.nls = predict(fit.nls, data.frame(v = v)),
        t.NA.pred.nls = ifelse(
            is.na(t.src),
            predict(fit.nls, data.frame(v = v)),
            t.src),
        t.pred.lm = exp(predict(fit.lm, data.frame(v = v)))) %>%
gather(source, t, 3:5) %>%
ggplot(aes(v, t, colour = source)) +
    geom_line() +
    geom_point(aes(v, t.src), colour = "black") +
    scale_colour_manual(values = c(
        "t.NA.pred.nls" = "black",
        "t.pred.nls" = "blue",
        "t.pred.lm" = "red"));

    enter image description here

    The black dots show the actual measurements. The blue curve is the power law nonlinear model fit, the red curve the linear fit after log-log transformation; the black curve corresponds to your original data, where v = NA values have been replaced with the nonlinear model fit estimates (so for v > 15 the black and blue curves overlap).