Nothing
predict.Mort2Dsmooth <-
function(object, newdata = NULL,
type = c("link", "response"),
se.fit = FALSE, ...){
## Input:
## object: a Mort2Dsmooth object
## newdata: an optional list in which
## to look for x and/or y
## variables with which to predict. If omitted,
## the fitted values are used.
## type: the type of prediction required.
## The default ("link") is
## on the scale of the linear predictors;
## the alternative "response" is on the
## scale of the response variable.
## se.fit: logical switch indicating if
## standard errors are required
## Output: a list with components
## fit: Predictions
## se.fit: Estimated standard errors
type <- match.arg(type)
## NO standard errors ## NO standard errors
if (!se.fit) {
## NO newdata ## NO newdata
if(is.null(newdata)) {
pred <- switch(type,
link = object$linear.predictors -
object$offset,
response = object$fitted.values)
}else{
## YES newdata ## YES newdata
## check newdata as a list
if(!is.list(newdata))
stop("newdata must be a data frame")
## check newdata names and variables
NEWdata <- list(x=object$x, y=object$y)
namesNEW <- names(NEWdata)
NEWdata[(namesnew <- names(newdata))] <- newdata
test <- noNms <- namesnew[!namesnew %in% namesNEW]
if(length(test) > 0){
stop("unknown names in newdata: ",
paste(noNms, collapse = ", "))
}
## original data
x <- object$x
y <- object$y
Z <- object$Z
offset <- object$offset
## new data
new.x <- sort(unique(c(x, newdata$x)))
new.y <- sort(unique(c(y, newdata$y)))
## new dimensions
new.m <- length(new.x)
new.n <- length(new.y)
## new weight matrix, where old and new coincide
new.W1 <- matrix(0, new.m, new.n)
new.W2 <- matrix(0, new.m, new.n)
whi.x <- which(new.x%in%x)
whi.y <- which(new.y%in%y)
new.W1[,whi.y] <- 10
new.W2[whi.x,] <- 1
whi <- (new.W1 - new.W2)==9
new.W <- matrix(0, new.m, new.n)
new.W[whi] <- 1
## new death and offset matrices
new.Z <- matrix(0, new.m, new.n)
new.Z[whi] <- Z
new.offset <- matrix(100, new.m, new.n)
new.offset[whi] <- offset
## new B-spline basis
## over x
if(min(new.x)>=min(x) & max(new.x)<=max(x)){
## only interpolation
new.ndx <- object$ndx[1]
xl <- min(new.x)
xr <- max(new.x)
xmax <- xr + 0.01 * (xr - xl)
xmin <- xl - 0.01 * (xr - xl)
new.Bx <- MortSmooth_bbase(x=new.x,
xl=xmin, xr=xmax,
ndx=new.ndx,
deg=object$deg[1])
}else{
## extrapolation + eventual interpolation
deg <- object$deg[1]
xl <- min(x)
xr <- max(x)
xmax <- xr + 0.01 * (xr - xl)
xmin <- xl - 0.01 * (xr - xl)
dx <- (xmax - xmin)/object$ndx[1]
xl1 <- min(new.x)
xr1 <- max(new.x)
minK0 <- xmin-deg:(deg+100)*dx
minK <- minK0[which(minK0<=xl1)[deg]]
maxK0 <- xmax+deg:(deg+100)*dx
maxK <- maxK0[which(maxK0>=xr1)[deg+1]]
knots <- seq(minK, maxK, by=dx)
PP <- outer(new.x, knots, MortSmooth_tpower, deg)
nn <- dim(PP)[2]
DD <-diff(diag(nn),diff=deg+1)/(gamma(deg+1)*dx^deg)
new.Bx <- (-1)^(deg + 1) * PP %*% t(DD)
}
nbx <- ncol(new.Bx)
## matplot(x, object$Bx, t="l", col=1,
## lty=1, xlim=range(new.x), lwd=2)
## matlines(new.x, new.Bx, t="l", col=2, lty=2, lwd=3)
## over y
if(min(new.y)>=min(y) & max(new.y)<=max(y)){
## only interpolation
new.ndx <- object$ndx[2]
yl <- min(new.y)
yr <- max(new.y)
ymax <- yr + 0.01 * (yr - yl)
ymin <- yl - 0.01 * (yr - yl)
new.By <- MortSmooth_bbase(x=new.y,
xl=ymin, xr=ymax,
ndx=new.ndx,
deg=object$deg[2])
}else{
deg <- object$deg[2]
yl <- min(y)
yr <- max(y)
ymax <- yr + 0.01 * (yr - yl)
ymin <- yl - 0.01 * (yr - yl)
dx <- (ymax - ymin)/object$ndx[2]
yl1 <- min(new.y)
yr1 <- max(new.y)
minK0 <- ymin-deg:(deg+100)*dx
minK <- minK0[which(minK0<=yl1)[deg]]
maxK0 <- ymax+deg:(deg+100)*dx
maxK <- maxK0[which(maxK0>=yr1)[deg+1]]
knots <- seq(minK, maxK, by=dx)
PP <- outer(new.y, knots, MortSmooth_tpower, deg)
nn <- dim(PP)[2]
DD <-diff(diag(nn),diff=deg+1)/(gamma(deg+1)*dx^deg)
new.By <- (-1)^(deg + 1) * PP %*% t(DD)
}
nby <- ncol(new.By)
## matplot(y, object$By, t="l", col=1,
## lty=1, xlim=range(new.y), lwd=2)
## matlines(new.y, new.By, t="l", col=2, lty=2, lwd=3)
## Row tensors of B-splines basis
Bx1 <- kronecker(matrix(1, ncol=nbx, nrow=1), new.Bx)
Bx2 <- kronecker(new.Bx, matrix(1, ncol=nbx,nrow=1))
RTBx <- Bx1*Bx2
By1 <- kronecker(matrix(1, ncol=nby, nrow=1), new.By)
By2 <- kronecker(new.By, matrix(1, ncol=nby, nrow=1))
RTBy <- By1*By2
## penalty stuff
Dx <- diff(diag(nbx), diff=object$pord[1])
Dy <- diff(diag(nby), diff=object$pord[2])
Px <- kronecker(diag(nby), t(Dx)%*%Dx)
Py <- kronecker(t(Dy)%*%Dy, diag(nbx))
## initialize
new.Z[is.na(new.Z)] <- 0
## 0) simple poisson-GLM with ages and years
## only for the interpolation cases
new.xx <- rep(new.x, new.n)
new.yy <- rep(new.y, each=new.m)
fit0 <- glm(round(c(new.Z)) ~ new.xx * new.yy +
offset(c(new.offset)),
family=poisson, weights=c(new.W))
etaGLM <- matrix(log(fit0$fitted) - c(new.offset),
new.m, new.n)
eta0 <- log((new.Z + 1)) - new.offset
eta0[new.W==0] <- etaGLM[new.W==0]
## simple ridge penalized regression
BBx <- solve(t(new.Bx)%*%new.Bx + diag(nbx) * 1e-6,
t(new.Bx))
BBy <- solve(t(new.By)%*%new.By + diag(nby) * 1e-6,
t(new.By))
a.init <- MortSmooth_BcoefB(BBx, BBy, eta0)
## estimating
new.fit <- Mort2Dsmooth_estimate(x=new.x,
y=new.y,
Z=new.Z,
offset=new.offset,
wei=new.W,
psi2=1,
Bx=new.Bx,
By=new.By,
nbx=nbx,nby=nby,
RTBx=RTBx, RTBy=RTBy,
lambdas=object$lambdas,
Px=Px,Py=Py,
a.init=a.init,
MON=FALSE,
TOL1=1e-6,
MAX.IT=50)
new.eta0 <- matrix(MortSmooth_BcoefB(new.Bx,
new.By,
new.fit$a),
new.m,new.n)
colnames(new.eta0) <- new.y
rownames(new.eta0) <- new.x
## matplot(x, log(Z/E), t="p", cex=0.5, pch=1,
## xlim=range(new.x), ylim=range(new.eta0))
## matlines(new.x, new.eta0, t="l", lty=2)
## matplot(y, t(log(Z/E)), t="p", cex=0.5, pch=1,
## xlim=range(new.y), ylim=range(new.eta0))
## matlines(new.y, t(new.eta0), t="l", lty=2)
## for(i in 1:new.m){
## bla <- new.eta0[i,]
## if(i%in%whi.x){
## bla1 <- log(Z[whi.x==i,]/E[whi.x==i,])
## }else{
## bla1 <- rep(NA,ncol(Z))
## }
## ran <- range(bla,bla1, na.rm=T)
## plot(new.y, t(bla), t="l", lty=2,
## xlim=range(new.y), ylim=ran,
## main=paste(new.x[i]))
## points(y, t(bla1))
## locator(1)
## }
## for(i in 1:new.n){
## bla <- new.eta0[,i]
## if(i%in%whi.y){
## bla1 <- log(Z[,whi.y==i]/E[,whi.y==i])
## }else{
## bla1 <- rep(NA,nrow(Z))
## }
## ran <- range(new.eta0)#bla,bla1, na.rm=T)
## plot(new.x, bla, t="l", lty=2,
## xlim=range(new.x), ylim=ran,
## main=paste(new.y[i]))
## points(x, t(bla1))
## locator(1)
## }
## in case either x or y is not provided in newdata
if(is.null(newdata$x)){
newdata$x <- x}
if(is.null(newdata$y)){
newdata$y <- y}
##
new.eta <- new.eta0[new.x%in%newdata$x,
new.y%in%newdata$y]
new.fitted0 <- exp(new.eta0 + new.offset)
new.fitted0[!whi] <- NA
new.fitted <- new.fitted0[new.x%in%newdata$x,
new.y%in%newdata$y]
pred <- switch(type,
link = new.eta,
response = new.fitted)
}
## YES standard errors ## YES standard errors
}else{
## NO newdata ## NO newdata
if(is.null(newdata)){
Bx <- object$Bx
By <- object$By
nbx <- ncol(Bx)
nby <- ncol(By)
## Row tensors of B-splines basis
Bx1 <- kronecker(matrix(1, ncol=nbx, nrow=1), Bx)
Bx2 <- kronecker(Bx, matrix(1, ncol=nbx,nrow=1))
RTBx <- Bx1*Bx2
By1 <- kronecker(matrix(1, ncol=nby, nrow=1), By)
By2 <- kronecker(By, matrix(1, ncol=nby, nrow=1))
RTBy <- By1*By2
## penalty stuff
Dx <- diff(diag(nbx), diff=object$pord[1])
Dy <- diff(diag(nby), diff=object$pord[2])
Px <- kronecker(diag(nby), t(Dx)%*%Dx)
Py <- kronecker(t(Dy)%*%Dy, diag(nbx))
P <- (object$lambdas[1]*Px) + (object$lambdas[2]*Py)
##
eta <- MortSmooth_BcoefB(Bx, By, object$coef)
mu <- exp(object$offset + eta)
W <- mu
WW <- object$W*W
BWB <- MortSmooth_BWB(RTBx, RTBy, nbx, nby, WW)
##
BWB.P1 <- solve(BWB + P)
se <- matrix(Mort2Dsmooth_se(RTBx=RTBx, RTBy=RTBy,
nbx=nbx, nby=nby,
BWB.P1=BWB.P1),
object$m, object$n,
dimnames=list(object$x, object$y))
## check usage of overdispersion
check.over <- eval(object$call$overdispersion)
if(is.null(check.over)){
check.over <- FALSE
}
if(check.over){
se.inflate <- sqrt(object$psi2)
se <- se * se.inflate
}
fit <- switch(type,
link = object$linear.predictors - object$offset,
response = object$fitted.values)
se.fit <- switch(type,
link = se,
response = object$fitted.values * (exp(se) -1))
}else{
## YES newdata ## YES newdata
## check newdata as a list
if(!is.list(newdata))
stop("newdata must be a data frame")
## check newdata names and variables
NEWdata <- list(x=object$x, y=object$y)
namesNEW <- names(NEWdata)
NEWdata[(namesnew <- names(newdata))] <- newdata
test <- noNms <- namesnew[!namesnew %in% namesNEW]
if(length(test) > 0){
stop("unknown names in newdata: ",
paste(noNms, collapse = ", "))
}
## original data
x <- object$x
y <- object$y
Z <- object$Z
offset <- object$offset
## new data
new.x <- sort(unique(c(x, newdata$x)))
new.y <- sort(unique(c(y, newdata$y)))
## new dimensions
new.m <- length(new.x)
new.n <- length(new.y)
## new weight matrix, where old and new coincide
new.W1 <- matrix(0, new.m, new.n)
new.W2 <- matrix(0, new.m, new.n)
whi.x <- which(new.x%in%x)
whi.y <- which(new.y%in%y)
new.W1[,whi.y] <- 10
new.W2[whi.x,] <- 1
whi <- (new.W1 - new.W2)==9
new.W <- matrix(0, new.m, new.n)
new.W[whi] <- 1
## new death and offset matrices
new.Z <- matrix(0, new.m, new.n)
new.Z[whi] <- Z
new.offset <- matrix(100, new.m, new.n)
new.offset[whi] <- offset
## new B-spline basis
## over x
if(min(new.x)>=min(x) & max(new.x)<=max(x)){
## only interpolation
new.ndx <- object$ndx[1]
xl <- min(new.x)
xr <- max(new.x)
xmax <- xr + 0.01 * (xr - xl)
xmin <- xl - 0.01 * (xr - xl)
new.Bx <- MortSmooth_bbase(x=new.x,
xl=xmin, xr=xmax,
ndx=new.ndx,
deg=object$deg[1])
}else{
## extrapolation + eventual interpolation
deg <- object$deg[1]
xl <- min(x)
xr <- max(x)
xmax <- xr + 0.01 * (xr - xl)
xmin <- xl - 0.01 * (xr - xl)
dx <- (xmax - xmin)/object$ndx[1]
xl1 <- min(new.x)
xr1 <- max(new.x)
minK0 <- xmin-deg:(deg+100)*dx
minK <- minK0[which(minK0<=xl1)[deg]]
maxK0 <- xmax+deg:(deg+100)*dx
maxK <- maxK0[which(maxK0>=xr1)[deg+1]]
knots <- seq(minK, maxK, by=dx)
PP <- outer(new.x, knots, MortSmooth_tpower, deg)
nn <- dim(PP)[2]
DD <-diff(diag(nn),diff=deg+1)/(gamma(deg+1)*dx^deg)
new.Bx <- (-1)^(deg + 1) * PP %*% t(DD)
}
nbx <- ncol(new.Bx)
## matplot(x, object$Bx, t="l", col=1,
## lty=1, xlim=range(new.x), lwd=2)
## matlines(new.x, new.Bx, t="l", col=2, lty=2, lwd=3)
## over y
if(min(new.y)>=min(y) & max(new.y)<=max(y)){
## only interpolation
new.ndx <- object$ndx[2]
yl <- min(new.y)
yr <- max(new.y)
ymax <- yr + 0.01 * (yr - yl)
ymin <- yl - 0.01 * (yr - yl)
new.By <- MortSmooth_bbase(x=new.y,
xl=ymin, xr=ymax,
ndx=new.ndx,
deg=object$deg[2])
}else{
deg <- object$deg[2]
yl <- min(y)
yr <- max(y)
ymax <- yr + 0.01 * (yr - yl)
ymin <- yl - 0.01 * (yr - yl)
dx <- (ymax - ymin)/object$ndx[2]
yl1 <- min(new.y)
yr1 <- max(new.y)
minK0 <- ymin-deg:(deg+100)*dx
minK <- minK0[which(minK0<=yl1)[deg]]
maxK0 <- ymax+deg:(deg+100)*dx
maxK <- maxK0[which(maxK0>=yr1)[deg+1]]
knots <- seq(minK, maxK, by=dx)
PP <- outer(new.y, knots, MortSmooth_tpower, deg)
nn <- dim(PP)[2]
DD <-diff(diag(nn),diff=deg+1)/(gamma(deg+1)*dx^deg)
new.By <- (-1)^(deg + 1) * PP %*% t(DD)
}
nby <- ncol(new.By)
## matplot(y, object$By, t="l", col=1,
## lty=1, xlim=range(new.y), lwd=2)
## matlines(new.y, new.By, t="l", col=2, lty=2, lwd=3)
## Row tensors of B-splines basis
Bx1 <- kronecker(matrix(1, ncol=nbx, nrow=1), new.Bx)
Bx2 <- kronecker(new.Bx, matrix(1, ncol=nbx,nrow=1))
RTBx <- Bx1*Bx2
By1 <- kronecker(matrix(1, ncol=nby, nrow=1), new.By)
By2 <- kronecker(new.By, matrix(1, ncol=nby, nrow=1))
RTBy <- By1*By2
## penalty stuff
Dx <- diff(diag(nbx), diff=object$pord[1])
Dy <- diff(diag(nby), diff=object$pord[2])
Px <- kronecker(diag(nby), t(Dx)%*%Dx)
Py <- kronecker(t(Dy)%*%Dy, diag(nbx))
## initialize
new.Z[is.na(new.Z)] <- 0
## 0) simple poisson-GLM with ages and years
## only for the interpolation cases
new.xx <- rep(new.x, new.n)
new.yy <- rep(new.y, each=new.m)
fit0 <- glm(round(c(new.Z)) ~ new.xx * new.yy +
offset(c(new.offset)),
family=poisson, weights=c(new.W))
etaGLM <- matrix(log(fit0$fitted) - c(new.offset),
new.m, new.n)
eta0 <- log((new.Z + 1)) - new.offset
eta0[new.W==0] <- etaGLM[new.W==0]
## simple ridge penalized regression
BBx <- solve(t(new.Bx)%*%new.Bx + diag(nbx) * 1e-6,
t(new.Bx))
BBy <- solve(t(new.By)%*%new.By + diag(nby) * 1e-6,
t(new.By))
a.init <- MortSmooth_BcoefB(BBx, BBy, eta0)
## estimating
new.fit <- Mort2Dsmooth_estimate(x=new.x,
y=new.y,
Z=new.Z,
offset=new.offset,
wei=new.W,
psi2=1,
Bx=new.Bx,
By=new.By,
nbx=nbx,nby=nby,
RTBx=RTBx, RTBy=RTBy,
lambdas=object$lambdas,
Px=Px,Py=Py,
a.init=a.init,
MON=FALSE,
TOL1=1e-6,
MAX.IT=50)
new.eta0 <- matrix(MortSmooth_BcoefB(new.Bx,
new.By,
new.fit$a),
new.m,new.n)
colnames(new.eta0) <- new.y
rownames(new.eta0) <- new.x
## in case either x or y is not provided in newdata
if(is.null(newdata$x)){
newdata$x <- x}
if(is.null(newdata$y)){
newdata$y <- y}
##
new.eta <- new.eta0[new.x%in%newdata$x,
new.y%in%newdata$y]
new.fitted0 <- exp(new.eta0 + new.offset)
new.fitted0[!whi] <- NA
new.fitted <- new.fitted0[new.x%in%newdata$x,
new.y%in%newdata$y]
## standard errors
BWB.P1 <- solve(new.fit$BWB + new.fit$P)
new.se0 <- matrix(Mort2Dsmooth_se(RTBx=RTBx,
RTBy=RTBy,
nbx=nbx, nby=nby,
BWB.P1=BWB.P1),
new.m, new.n,
dimnames=list(new.x, new.y))
new.se <- new.se0[new.x%in%newdata$x,
new.y%in%newdata$y]
## check usage of overdispersion
check.over <- eval(object$call$overdispersion)
if(is.null(check.over)){
check.over <- FALSE
}
if(check.over){
se.inflate <- sqrt(object$psi2)
new.se <- new.se * se.inflate
}
fit <- switch(type,
link = new.eta,
response = new.fitted)
se.fit <- switch(type,
link = new.se,
response =new.fitted*(exp(new.se)-1))
}
pred <- list(fit=fit, se.fit=se.fit)
}
return(pred)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.