Next: Examples Up: ASSIST: A Suite of Previous: The snm Function

# Vector Spline Models

In many applications two or more dependent variables are observed at several values of independent variables such as at multiple time points. Often observations from different variables are contemporaneously correlated. Observations from the same variable may also be correlated. The statistical problems are (i) to estimate functions that model their dependences on the independent variables and (ii) to investigate relationships between these functions. Wang et al. (2000) proved that the joint estimates have smaller posterior variances than those of function-by-function estimates and are therefore more efficient. In this section we show how to use ssr and nnr to fit functions jointly.

Consider the following model

 (48)

where the th response of the th variable is generated as the th function evaluated at the design point plus a random error . We assume that has the model space ().

Denote , , , , , , and . Assume that , where depends on some parameters .

We estimate all functions 's jointly as the minimizer of the following penalized weighted least squares

 (49)

where is the orthogonal projection operator of onto in . We use the GML method to estimate the variance-covariance parameters and the smoothing parameters 's as the minimizers of ().

We now show how to trick ssr to fit model (). For the simplicity of notations, we consider the case of . Situations with can be fitted similarly. Note that the domains of may be different for different . For most applications they are the same, which is assumed in the remaining of this section. Denote the common domain as . Rewrite as , which is now considered as a function of both and variables on the tensor product domain . Then we can represent the original functions as

 (50)

Let
 (51)

be the model space for , where , and 's are RKHS's with RK for . Then it is easy to check that
 (52)

where , for , and for . are RKHS's with RKs for and for . The model space () is similar to that of a SS ANOVA model. Thus we can use the function ssr to fit functions 's jointly.

We may reparametrize as

 (53) (54)

() and () are SS ANOVA decompositions of with the set-to-zero and sum-to-zero side conditions respectively. This kind of ANOVA decomposition can be carried out for general .

For (), let and . Let in () be the model space of . Then

 (55)

where , for , and for . are RKHS's with RKs for and for .

For (), denote and . Let in () be the model space of . Then

 (56)

where , for , and for . are RKHS's with RKs for and for .

Often we are interested in possible relationships, if any, between and . For example, one may want to check if they are equal or parallel. Let be a probability measure on . Consider the following SS ANOVA decomposition

 (57)

where

We have . is the overall mean, and are the main effects and is the interaction. (), equivalent of (), will make some hypotheses more transparent. It is easy to check that the following hypotheses are equivalent:

Furthermore, if and ,

Therefore the hypothesis that and are equal is equivalent to the effect . The hypothesis that and are parallel is equivalent to the hypothesis that the interaction . The hypothesis that the integral of equals to that of is equivalent to the hypothesis that the main effect . The hypothesis that the summation of and is a constant is equivalent to the hypothesis that the main effect . The hypothesis that there exists a linear relationship between the functions and is equivalent to the hypothesis that the interaction is multiplicative. Thus, for these simple hypotheses we can fit the SS ANOVA model () and perform tests on the corresponding components.

For illustration, we generate a data set from the following model

where random errors 's follow bivariate normal with , and . The variance-covariance structure can be specified with the combination of the weights and correlation options. We will fit with an arbitrary pairwise variance-covariance structure.

Suppose we want to use cubic splines to model and in (), and in (), and and in (). Then , , , , and and are the RK of a cubic spline. In the following we first fit and using marginal data as bisp.fit1 and bisp.fit2 respectively. Then we fit jointly using the formulation () as bisp.fit3 the formulation () as bisp.fit4, and the formulation (), or equivalently the formulation (), as bisp.fit5.

> options(contrasts=rep("contr.treatment", 2))
> n <- 100
> s1 <- .5
> s2 <- 1
> r <- .5
> A <- diag(c(s1,s2))%*%matrix(c(sqrt(1-r**2),0,r,1),2,2)
> e <- NULL
> for (i in 1:n) e <- c(e,A%*%rnorm(2))
> t <- 1:n/n
> y1 <- sin(2*pi*t) + e[seq(1,2*n,by=2)]
> y2 <- sin(2*pi*t) + 2*t + e[seq(2,2*n,by=2)]
> bisp.dat <- data.frame(y=c(y1,y2),t=rep(t,2),id=rep(c(0,1),rep(n,2)),
pair=rep(1:n,2))

# fit separately
> bisp.fit1 <- ssr(y~I(t-.5),rk=cubic(t),spar="m",
data=bisp.dat[bisp.dat$id==0,]) > p.bisp.fit1 <- predict(bisp.fit1) > bisp.fit2 <- ssr(y~I(t-.5),rk=cubic(t),spar="m", data=bisp.dat[bisp.dat$id==1,])
> p.bisp.fit2 <- predict(bisp.fit2)

# fit jointly
> bisp.fit3 <- ssr(y~id*I(t-.5), rk=list(rk.prod(cubic(t),kron(id==0)),
rk.prod(cubic(t),kron(id==1))), spar="m",
weights=varIdent(form=~1|id),
cor=corSymm(form=~1|pair), data=bisp.dat)
> bisp.fit3
...
GML estimate(s) of smoothing parameter(s) : 0.2793703 0.3788567
Equivalent Degrees of Freedom (DF):  11.84544
Estimate of sigma:  0.4616973
Correlation structure of class corSymm representing
Correlation:
1
2 0.523
Variance function structure of class varIdent representing
0        1
1 2.270982

> p.bisp.fit3.1 <- predict(bisp.fit3,newdata=bisp.dat[bisp.dat$id==0,], terms=c(1,0,1,0,1,0)) > p.bisp.fit3.2 <- predict(bisp.fit3,newdata=bisp.dat[bisp.dat$id==1,],
terms=c(1,1,1,1,0,1))
> p.bisp.fit3.1$pstd <- p.bisp.fit3.1$pstd*sqrt((2*n-bisp.fit3$df) /(2*n-bisp.fit3$df-1))
> p.bisp.fit3.2$pstd <- p.bisp.fit3.2$pstd*sqrt((2*n-bisp.fit3$df) /(2*n-bisp.fit3$df-1))

> bisp.fit4 <- ssr(y~id*I(t-.5), rk=list(cubic(t),
rk.prod(cubic(t),kron(id==1))),
spar="m", weights=varIdent(form=~1|id),
cor=corSymm(form=~1|pair), data=bisp.dat)
> p.bisp.fit4.1 <- predict(bisp.fit4,newdata=bisp.dat[bisp.dat$id==0,], terms=c(1,0,1,0,1,0)) > p.bisp.fit4.2 <- predict(bisp.fit4,newdata=bisp.dat[bisp.dat$id==1,],
terms=c(1,1,1,1,1,1))
> p.bisp.fit4.3 <- predict(bisp.fit4,newdata=bisp.dat[bisp.dat$id==1,], terms=c(0,1,0,1,0,1)) > p.bisp.fit4.4 <- predict(bisp.fit4,newdata=bisp.dat[bisp.dat$id==1,],
terms=c(0,0,0,1,0,1))
> p.bisp.fit4.1$pstd <- p.bisp.fit4.1$pstd*sqrt((2*n-bisp.fit4$df) /(2*n-bisp.fit4$df-1))
> p.bisp.fit4.2$pstd <- p.bisp.fit4.2$pstd*sqrt((2*n-bisp.fit4$df) /(2*n-bisp.fit4$df-1))
> p.bisp.fit4.3$pstd <- p.bisp.fit4.3$pstd*sqrt((2*n-bisp.fit4$df) /(2*n-bisp.fit4$df-1))
> p.bisp.fit4.4$pstd <- p.bisp.fit4.4$pstd*sqrt((2*n-bisp.fit4$df) /(2*n-bisp.fit4$df-1))
> bisp.fit5 <- ssr(y~id*I(t-.5),
rk=list(cubic(t),rk.prod(cubic(t),kron((id==0)-(id==1)))),
spar="m", weights=varIdent(form=~1|id),
cor=corSymm(form=~1|pair), data=bisp.dat)
> p.bisp.fit5.1 <- predict(bisp.fit5,newdata=bisp.dat[bisp.dat$id==0,], terms=c(1,0,1,0,1,1)) > p.bisp.fit5.2 <- predict(bisp.fit5,newdata=bisp.dat[bisp.dat$id==1,],
terms=c(1,1,1,1,1,1))
> p.bisp.fit5.3 <- predict(bisp.fit5,newdata=bisp.dat[bisp.dat$id==1,], terms=c(0,1,0,1,0,1)) > p.bisp.fit5.4 <- predict(bisp.fit5,newdata=bisp.dat[bisp.dat$id==1,],
terms=c(0,0,0,1,0,1))
> p.bisp.fit5.1$pstd <- p.bisp.fit5.1$pstd*sqrt((2*n-bisp.fit5$df) /(2*n-bisp.fit5$df-1))
> p.bisp.fit5.2$pstd <- p.bisp.fit5.2$pstd*sqrt((2*n-bisp.fit5$df) /(2*n-bisp.fit5$df-1))
> p.bisp.fit5.3$pstd <- p.bisp.fit5.3$pstd*sqrt((2*n-bisp.fit5$df) /(2*n-bisp.fit5$df-1))
> p.bisp.fit5.4$pstd <- p.bisp.fit5.4$pstd*sqrt((2*n-bisp.fit5$df) /(2*n-bisp.fit5$df-1))


For each fit, we calculate the estimated functions and posterior variances. For fits bisp.fit4 and bisp.fit5, we also calculate the estimates and posterior variances of and in (). They are used to check if and if they are parallel respectively. Note that we inflate the posterior variances by one more degree of freedom used for estimating the correlation parameter. These estimates are shown in Figures , and . Even though not obvious in Figure , the Bayesian confidence intervals of the joint estimates are uniformly narrower than those of the function-by-function estimates. Obviously from Figures and that these two functions are not equal, nor parallel.

Next: Examples Up: ASSIST: A Suite of Previous: The snm Function
Yuedong Wang 2004-05-19