Fitting two sets of values by changing x-axis

Question

I have two datasets. The first one consist of measurements with a known date. The second dataset consists of measurements along a distance axis. The second dataset is sequentially deposited (for example an ice-core) with an unknown growth rate. Values in the two datasets are expected to match within a certain random error. An example:

curve <- data.frame(date = as.Date(c('01/09/11', '01/10/11', '01/11/11', '01/12/11', '01/01/12', '01/02/12', '01/03/12', '01/04/12', '01/05/12', '01/06/12', '01/07/12', '01/08/12', '01/09/12'), format = "%d/%m/%y"), value = c(3.33, 3.80, 4.19, 4.43, 4.56, 4.76, 4.71, 4.96, 4.72, 4.66, 4.35, 4.12, 3.88))

points <- data.frame(distance = 1:10, value = c(3.5, 4.2, 4.5, 4.4, 4.7, 4.8, 5.1, 4.9, 4.1, 3.7))

library(ggplot2)
library(scales)

a <- ggplot(curve, aes(x = date, y = value)) + geom_line() + scale_x_date(name = "Month", breaks = date_breaks("months"), labels = date_format("%b"))
b <- ggplot(points, aes(x = distance, y = value)) + geom_point()

grid.arrange(a,b,ncol=1)

I call the dataset with known date "curve" because the number of measurements is much larger than for the dataset with a distance axis ("point" dataset). In reality both of them consist of measurement points, but I suppose I could approximate the curve dataset to an actual curve by using loess or splines. I would like to fit the "point" dataset to the "curve" dataset by changing the x-axis (or growth rate) for the "point" dataset with a constraint that the growth rate > 0 (distance between "point"s). I assume that I would be making a growth model for the distance axis.

What would be the best method to do this kind of fitting? Is there an R package developed for this kind of problem?

Answer 1

@Glen_b suggested using dynamic time warping. There is an R package for DTW calculations. It works for the example, except for that I have not managed to implement the constraint that "growth rate" or rather minimum step in this context should be > 0. On the other hand implementing such a constraint for the example data would not make much sense because of few reference index points. So here is how it goes:

library(dtw)
res <- dtw(points$value, curve$value, keep = T, step=typeIIIc)

library(plyr)
names <- data.frame(date = curve$date, 
              index = seq_along(curve$date))

points$date <- as.Date(mapvalues(res$index2, from = names$index, 
        to = as.character(names$date)))

library(ggplot2)
library(scales)

ggplot() + geom_line(data = curve, aes(x = date, y = value)) + 
 geom_text(data = points, aes(x=date, y=value, label = distance), 
 color = "red") + scale_x_date(name = "Month", 
 breaks = date_breaks("months"), labels = date_format("%b"))

Red numbers are the fitted points. Line represents the "curve". Let´s try this on more reference index points:

new.date <- seq.Date(from = curve$date[1], 
      to = curve$date[length(curve$date)], by = "day")

curve <- merge(curve, data.frame(date = new.date, 
     index = seq_along(new.date)), by = "date", all = F, sort = F)
lo.c <- loess(value ~ index, curve, span = 0.4)

new.curve <- data.frame(date = new.date, 
       index = seq_along(new.date), 
       value = predict(lo.c, seq_along(new.date)))

## New.curve
ggplot() + geom_line(data = new.curve, aes(x = date, y = value)) + 
    geom_point(data = curve, aes(x = date, y = value), 
    color = "red") + scale_x_date(name = "Month", 
    breaks = date_breaks("months"), labels = date_format("%b")) 

Red points are the old reference index points. Line represents loess curve fitted on the "curve" values.

I encounter problems when trying to define the step for dtw:

res.new <- dtw(points$value, new.curve$value, keep = T, 
       step=typeIIIc)

#Error in dtw(points$value, new.curve$value, keep = T, 
  step = typeIIIc) : 
      No warping paths exists that is allowed by constraints

If I remove the step, fitting works, but now the function fits query indices on reference indices n = length(new.date).

res.new <- dtw(points$value, new.curve$value, keep = T)
plot(res.new, "twoway") 

I want only one time point for each query. Not an average over a longer period, so I am having difficulties to make this method work for the problem. This is no doubt my fault, because dtw package documentation is not very easy to understand. Any suggestions how to proceed?