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)

enter image description here

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?

  • 3
    It sounds a bit like you're asking for something akin to dynamic time warping. That garners a few hits here as a search term, though there are a number of other ways to try to estimate some continuous nonlinear relationship between your two variables. – Glen_b Dec 3 '13 at 0:21
  • I would have guessed that you could fairly easily use nls or one of the other optimizers in R to solve for the offset to the points that would minimize the (value-value)-differences between the curve and the point-set. It would require that we know that the intervals are on the same scale. – DWin Dec 3 '13 at 0:41
  • Dynamic time wrapping sounds promising. Thanks for pointing that out. I will look into it. There are papers on this topic, for instance: link.springer.com/article/10.1007%2Fs11004-011-9352-7. They use an approach which is probably connected to @Dwin´s suggestion. I am not a mathematician, so understanding what they did takes time... – Mikko Dec 3 '13 at 8:27

@Glen_b suggested using dynamic time wrapping. 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"))

enter image description here

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 costraints

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

enter image description here

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?

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.