Concurrent effects:

\( y_i(t) = f(x_i(t)) + e_i(t) \)

library(refund)

## Loading required package: fda
## Loading required package: splines
## Loading required package: Matrix
## 
## Attaching package: 'fda'
## 
## The following object is masked from 'package:graphics':
## 
##     matplot

Generate data

n <- 50
g <- 100
t_vec <- seq(.1, .9, l=g)
sd_error <- .3

set.seed(20140620)
x <- t(replicate(n, rnorm(1) + 
                   dbeta(t_vec, runif(1, 1, 10), runif(1,1,10))))
# substract mean x-trajectory for centered covariate: 
# (not strictly necessary for real data, only here to avoid level shifts in 
# between estimated functions and functions used to generate the data):
x <- sweep(x, 2, colMeans(x))

data <- data.frame(id=1:n)
data$x <- x

Generate \( y(t) \)

\( f(x(t)) = x(t)\beta(t) \) with \( \beta(t) \equiv 1 \):

data$yconst <- x + sd_error * matrix(rnorm(n * g), nrow=n)

\( f(x(t)) = x(t) \log(t) \); i.e. \( \beta(t) = \log(t) \):

data$ylin <- t(t(x) * log(t_vec)) + sd_error * matrix(rnorm(n * g), nrow=n)

\( f(x(t)) = (x(t)/4)^3 \)

data$ysmoo <- (x/4)^3 + sd_error * matrix(rnorm(n * g), nrow=n)

\( f(x(t)) = \sin(t x) \)

data$ysmoosmoo <- t(sin(t_vec * t(x))) + 
    sd_error * matrix(rnorm(n * g), nrow=n)

that's what it looks like:

layout(t(matrix(1:6, 3, 2)))
matplot(t_vec, t(data$x[1:10,]), type="l", lty=1, col=rgb(0,0,0,.1))
matplot(t_vec, t(data$yconst[1:10,]), type="l", lty=1, col=rgb(0,0,0,.1))
matplot(t_vec, t(data$ylin[1:10,]), type="l", lty=1, col=rgb(0,0,0,.1))
matplot(t_vec, t(data$ysmoo[1:10,]), type="l", lty=1, col=rgb(0,0,0,.1))
matplot(t_vec, t(data$ysmoosmoo[1:10,]), type="l", lty=1, col=rgb(0,0,0,.1))

Fit models & check estimates

~ c(x) : linear effect of x constant over index of y

mconst <- pffr(yconst ~ c(x), data=data, yind=t_vec) 
summary(mconst)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## yconst ~ c(x)
## 
## Constant coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.00173    0.00434    -0.4     0.69    
## x            0.99695    0.00345   289.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Smooth terms & functional coefficients:
##                      edf Ref.df F p-value
## Intercept(t_vec) 0.00197     19 0    0.75
## 
## R-sq.(adj) =  0.944   Deviance explained = 94.4%
## REML score = 1193.2  Scale est. = 0.094052  n = 5000 (50 x 100)

coefficient of x is approx. 1

~ x: linear effect of x varying over index of y:

mlin <- pffr(ylin ~ x, data=data, yind=t_vec)

layout(1)
plot(mlin, select=2) #plot beta(t)
curve(log(x), lty=2, col="red", add=TRUE)

…looks good

~ c(s(x)): smooth effect of x constant over index of y

msmoo <- pffr(ysmoo ~ c(s(x)), data=data, yind=t_vec) 

plot(msmoo, select=2)
curve((x/4)^3, lty=2, col="red", add=TRUE)

…looks good as well

~ s(x): smooth effect of x varying over index of y:

msmoosmoo <- pffr(ysmoosmoo ~ s(x), data=data, yind=t_vec) 
x_grid <- seq(min(x), max(x), l=40)
t_grid <- seq(min(t_vec), max(t_vec), l=40)
true_smoosmoo <- {
    grid <- expand.grid(x=x_grid, t=t_grid)
    f <- with(grid, sin(t*x))
    matrix(f, 40, 40)
}

layout(t(matrix(1:4, 2, 2))); 
plot(msmoosmoo, select=2, se=FALSE, rug=FALSE, scheme=2)
image(x=x_grid, y=t_grid, z=true_smoosmoo, main="true effect")
contour(x=x_grid, y=t_grid, z=true_smoosmoo, add=TRUE)
plot(msmoosmoo, select=2, scheme=1, 
     ticktype="detailed", zlim=range(true_smoosmoo))

## Warning: surface extends beyond the box

persp(x=x_grid, y=t_grid, z=true_smoosmoo, 
      ticktype="detailed", zlim=range(true_smoosmoo), theta = 30, phi = 30)

not too bad, either.