R/spmDesign.r In SemiPar: Semiparametic Regression

```########## S-function: spmDesign ##########

# For creating the design matrices for
# a semiparametric model.

# Last changed: 22 SEP 2005

spmDesign <- function(spm.info)
{
# Extract required information

degree.pen <- spm.info\$pen\$degree
degree.krige <- spm.info\$krige\$degree
knots.pen <- spm.info\$pen\$knots
basis.pen <- spm.info\$pen\$basis
x.pen <- spm.info\$pen\$x

# Initialise matrices and counters

X <- NULL
Z <- NULL
block.inds <- NULL
re.block.inds <- NULL
trans.mat <- NULL
stt.pos <- 1
re.stt.pos <- 1

# Obtain X-matrix components corresponding
# to the "block.inds", starting with the
# intercept.

block.inds <- 1
stt.pos <- stt.pos + 1

intercept.only <- (is.null(spm.info\$lin)&is.null(spm.info\$pen)
&is.null(spm.info\$krige))

if (!is.null(spm.info\$lin))
{
for (j in 1:ncol(as.matrix(spm.info\$lin\$x)))
block.inds <- c(block.inds,list(j+stt.pos-1))

stt.pos <- stt.pos + ncol(as.matrix(spm.info\$lin\$x))
}

if (!is.null(spm.info\$pen))
{

ncol.X.pen <- degree.pen
if(basis.pen=="tps")
{
ncol.X.pen <- (ncol.X.pen-1)/2
if (any(floor(ncol.X.pen)!=ncol.X.pen))
stop("Only odd degree thin plate splines supported.")
}

for (j in 1:ncol(as.matrix(spm.info\$pen\$x)))
{
block.inds <- c(block.inds,list(stt.pos:(stt.pos+ncol.X.pen[j]-1)))
stt.pos <- stt.pos + ncol.X.pen[j]
}

}

if (!is.null(spm.info\$krige))
{
ncol.X.krige <- (degree.krige+2)*(degree.krige+4)/8 - 1

block.inds <- c(block.inds,list(stt.pos:(stt.pos+ncol.X.krige-1)))
stt.pos <- stt.pos + ncol.X.krige
}

# Now obtain the X and Z matrices, and update
# "block.inds", "re.block.inds" and "trans.mat" sequentially.

comp.num <- 2
if (!is.null(spm.info\$lin))
{
# Update the X matrix

X <- cbind(X,spm.info\$lin\$x)

comp.num <- comp.num + ncol(X)
}

if (!is.null(spm.info\$pen))
{

for (j in 1:ncol(as.matrix(x.pen)))
{
# Update the X matrix.

for (ideg in 1:ncol.X.pen[j])
X <- cbind(X,as.matrix(x.pen)[,j]^ideg)

# Update the Z matrix.

if (basis.pen=="tps")
{
new.cols <- outer(as.matrix(x.pen)[,j],knots.pen[[j]],"-")
new.cols <- abs(new.cols)^degree.pen[j]

sqrt.Omega <- matrix.sqrt(abs(outer(knots.pen[[j]],
knots.pen[[j]],"-"))^degree.pen[j])

new.cols <- t(solve(sqrt.Omega,t(new.cols)))
trans.mat <- c(trans.mat,list(sqrt.Omega))
}

if (basis.pen=="trunc.poly")
{
new.cols <- outer(as.matrix(x.pen)[,j],knots.pen[[j]],"-")
new.cols <- (new.cols*(new.cols>0))^degree.pen[j]
trans.mat <- c(trans.mat,list(NULL))
}

Z <- cbind(Z,new.cols)

re.block.inds <- c(re.block.inds,list(re.stt.pos:ncol(Z)))

end.pos <- ncol(Z) + stt.pos - re.stt.pos
block.inds[[comp.num]] <- c(block.inds[[comp.num]],stt.pos:end.pos)

re.stt.pos <- ncol(Z) + 1
stt.pos <-  end.pos + 1
comp.num <- comp.num + 1
}
}

if (!is.null(spm.info\$krige))
{
# Now add on the kriging contribution, starting with
# obtaining the "Omega matrix" of intra-knot covariances,
# but first work out range parameter.

knots.krige <- spm.info\$krige\$knots

# Update the X matrix

x1.v <- spm.info\$krige\$x[,1]
x2.v <- spm.info\$krige\$x[,2]

m.krige <- (degree.krige/2) + 1

for (im in 2:m.krige)
{
X <- cbind(X,x1.v^(im-1))
if (im>2)
for (ipow in 2:(im-1))
X <- cbind(X,(x1.v^(im-ipow))*(x2.v^(ipow-1)))
X <- cbind(X,x2.v^(im-1))
}

# Update the Z matrix

num.knots.krige <- nrow(knots.krige)

dist.mat <- matrix(0,num.knots.krige,num.knots.krige)
dist.mat[lower.tri(dist.mat)] <- dist(as.matrix(knots.krige))
dist.mat <- dist.mat + t(dist.mat)

Omega <- tps.cov(dist.mat,m=m.krige,d=2)

# Obtain preliminary Z matrix of knot to data covariances,
# for the kriging component.

x.knot.diffs.1 <- outer(spm.info\$krige\$x[,1],knots.krige[,1],"-")
x.knot.diffs.2 <- outer(spm.info\$krige\$x[,2],knots.krige[,2],"-")
x.knot.dists <- sqrt(x.knot.diffs.1^2+x.knot.diffs.2^2)

prelim.Z <- tps.cov(x.knot.dists,m=m.krige,d=2)

# Transform to mixed model canonical form
# (covariance matrix of random effects is multiple
# of identity)

sqrt.Omega <- matrix.sqrt(Omega)

trans.mat <- c(trans.mat,list(sqrt.Omega))

# Combine to form final Z matrix

Z <- cbind(Z,t(solve(sqrt.Omega,t(prelim.Z))))

re.block.inds <- c(re.block.inds,list(re.stt.pos:ncol(Z)))

end.pos <- ncol(Z) + stt.pos - re.stt.pos
block.inds[[comp.num]] <- c(block.inds[[comp.num]],stt.pos:end.pos)

re.stt.pos <- ncol(Z) + 1
stt.pos <- end.pos + 1
comp.num <- comp.num + 1
}

# Add on column of 1's for the intercept term

if (intercept.only)
X <- as.matrix(rep(1,length(spm.info\$y)))

if (!intercept.only)
X <- cbind(rep(1,nrow(X)),X)

# Save spline component of Z separately.

Z.spline <- Z

# If random term present then add on Kronecker-type
# intercept structure to Z matrix

if (!is.null(spm.info\$random))
{
group.inds <- spm.info\$random\$group.inds

num.groups <- length(group.inds)

for (j in 1:num.groups)

# Append to end of Z matrix

re.block.inds <- c(re.block.inds,list(re.stt.pos:ncol(Z)))

end.pos <- ncol(Z) + stt.pos - re.stt.pos

block.inds <- c(block.inds,list(stt.pos:end.pos))

trans.mat <- c(trans.mat,list(NULL))
}

if (is.null(spm.info\$krige)) sqrt.Omega.krige <- NULL

# Return spm() fit object

return(list(spm.info=spm.info,X=X,Z=Z,Z.spline=Z.spline,
trans.mat=trans.mat,block.inds=block.inds,
re.block.inds=re.block.inds))

}

########## End of spmDesign ##########
```

Try the SemiPar package in your browser

Any scripts or data that you put into this service are public.

SemiPar documentation built on May 2, 2019, 5:42 a.m.