In larslau/pavpop: Phase and amplitude varying population pattern models

knitr::opts_chunk$set(warning = FALSE, comment = "#>", fig.height = 5, fig.width = 7)
library(pavpop)

Generate data

# Number of samples
n <- 20 
# Number of observation points
m <- 100 

# Observation points
t <- seq(0, 1, length = m + 2)[2:(m + 1)] 

# Common basis function (both mean and amplitude variation)
kts <- seq(0, 1, length = 10)
basis_fct <- make_basis_fct(kts = kts, intercept = TRUE, 
                            control = list(boundary = c(-0.5, 1.5)))
amp_fct <- make_basis_fct(type = 'intercept')
df <- attr(basis_fct, 'df')

# Generate true mean weights 
beta_t <- rexp(df, 0.5)

# Generate random intercepts
b_t <- rnorm(n, sd = 2)

# Generate warping function and random parameters
sigma <- 0.1
warp_fct <- make_warp_fct(type = 'linear')
warp_cov_t <- matrix(c(2, 0.5, 0.5, 1), 2, 2)
w_t <- replicate(n, (t(chol(warp_cov_t)) %*% rnorm(2, sd = sigma))[, 1])

# Generate data
y <- lapply(1:n, function(i) {as.numeric(basis_fct(warp_fct(w_t[, i], t)) %*% beta_t 
                                         + amp_fct(t) %*% b_t[i] 
                                         + rnorm(m, sd = sigma))}) 
t <- lapply(1:n, function(x) t)

# Plot observations
plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', 
     xlab = 't', ylab = 'y(t)')
legend(0.7, range(y)[2], legend = expression(theta(t)), lty = 2, lwd = 2)
for (i in 1:n) lines(t[[i]], y[[i]], col = rainbow(n)[i])
lines(t[[1]], basis_fct(t[[1]]) %*% beta_t , lwd = 2, lty = 2)

pavpop estimation

We now set up pavpop to estimate in the sitar model

amp_cov <- make_cov_fct(id_cov, noise = FALSE)
warp_cov <- make_cov_fct(unstr_cov, param = c(1, 1, 0), noise = FALSE)

res <- pavpop(y, t, basis_fct, warp_fct, amp_cov, warp_cov, amp_fct, iter = c(10, 10))

# Plot aligned samples
t_p <- seq(-0.2, 1.2, length = m + 2)
plot(0, 0, xlim = c(-0.2, 1.2), ylim = range(y), type = 'n', 
     xlab = 'warped t', ylab = 'y(t)')
legend(1, range(y)[2], legend = c(expression(theta(t)), expression(hat(theta)(t))), 
       lty = 2, lwd = 2, col = c('black', 'red'))
for (i in 1:n) lines(warp_fct(res$w[, i], t[[i]]), y[[i]], col = rainbow(n)[i])
lines(t_p, basis_fct(t_p) %*% beta_t , lwd = 2, lty = 2)
lines(t_p, basis_fct(t_p) %*% res$c, lwd = 2, lty = 2, col = 'red')

plot(as.numeric(w_t), as.numeric(res$w), xlab = 'true warp variables', 
     ylab = 'predicted warp variables', pch = 19, col = c('green', 'blue'), asp = 1)
abline(0, 1, lty = 2)

Are the parameter estimates okay?

# Noise standard deviation: 0.1
res$sigma
# Amplitude standard deviations: 2
res$sigma * sqrt(res$amp_cov_par)

# Warp covariance

# True
sigma^2 * warp_cov_t
# Observed
cov(t(w_t))
# estimated
res$sigma^2 * warp_cov(1:2, res$warp_cov_par)

Traditional sitar estimation

library(sitar)

dat <- data.frame(x = unlist(t), y = unlist(y), id = rep(1:n, each = m))

res_sitar <- sitar(x = x, y = y, id = id, data = dat, df = 14)


plot(res_sitar, opt = 'd', xlim = c(-0.2, 1.2), ylim = range(y), lwd = 2, lty = 2, 
     col = 'red', xlab = 'warped t', ylab = 'y(t)')
lines(t_p, basis_fct(t_p) %*% beta_t , lwd = 2, lty = 2)
legend(1, range(y)[2], legend = c(expression(theta(t)), expression(hat(theta)(t))), 
       lty = 2, lwd = 2, col = c('black', 'red'))


# Plot warped curves
abc <- random.effects(res_sitar)
abc$a <- 0
warped <- with(dat, xyadj(x, y, id, res_sitar, abc))
for (i in 1:n) 
  with(warped, lines(x[1:m + (i - 1) * m], y[1:m + (i - 1) * m], col = rainbow(n)[i]))

plot(as.numeric(w_t), rbind(-abc$b * exp(abc$c), exp(abc$c) - 1), 
     xlab = 'true warp variables', ylab = 'predicted warp variables', 
     pch = 19, col = c('green', 'blue'), asp = 1)
abline(0, 1, lty = 2)