.egamFit <- function(x, y, sp, Eb, UrS,
weights, start, offset, U1, Mp, family,
control, null.coef, needVb) {
## Routine for fitting GAMs beyond exponential family.
## Inputs as gam.fit3 except that family is of class "extended.family", while
## sp contains the vector of extended family parameters, followed by the log smoothing parameters,
scoreType <- "REML"
deriv <- 0
theta <- family$getTheta()
## penalized <- if (length(UrS)>0) TRUE else FALSE
x <- as.matrix(x)
nSp <- length(sp)
rank.tol <- .Machine$double.eps*100 ## tolerance to use for rank deficiency
q <- ncol(x)
n <- nobs <- nrow(x)
xnames <- dimnames(x)[[2]]
ynames <- if (is.matrix(y)) rownames(y) else names(y)
## Now a stable re-parameterization is needed....
if (length(UrS)) {
rp <- gam.reparam(UrS, sp, deriv)
T <- diag(q)
T[1:ncol(rp$Qs),1:ncol(rp$Qs)] <- rp$Qs
T <- U1%*%T ## new params b'=T'b old params
null.coef <- t(T)%*%null.coef
# Start is a list of vectors, each is a different possible initialization
for(jj in 1:length(start)){ start[[jj]] <- t(T)%*%start[[jj]] }
## form x%*%T in parallel
x <- .Call(.C_qgam_pmmult2,x,T,0,0,control$nthreads)
rS <- list()
for (i in 1:length(UrS)) {
rS[[i]] <- rbind(rp$rS[[i]],matrix(0,Mp,ncol(rp$rS[[i]])))
} ## square roots of penalty matrices in current parameterization
Eb <- Eb%*%T ## balanced penalty matrix
rows.E <- q-Mp
Sr <- cbind(rp$E,matrix(0,nrow(rp$E),Mp))
St <- rbind(cbind(rp$S,matrix(0,nrow(rp$S),Mp)),matrix(0,Mp,q))
} else {
T <- diag(q);
St <- matrix(0,q,q)
rSncol <- rows.E <- Eb <- Sr <- 0
rS <- list(0)
rp <- list(det=0,det1 = 0,det2 = 0,fixed.penalty=FALSE)
}
## re-parameterization complete. Initialization....
nvars <- ncol(x)
if (nvars==0) stop("emtpy models not available")
if (is.null(weights)) weights <- rep.int(1, nobs)
if (is.null(offset)) offset <- rep.int(0, nobs)
linkinv <- family$linkinv
valideta <- family$valideta
validmu <- family$validmu
dev.resids <- family$dev.resids
## need an initial `null deviance' to test for initial divergence...
## if (!is.null(start)) null.coef <- start - can be on edge of feasible - not good
null.eta <- as.numeric(x%*%null.coef + as.numeric(offset))
## call the families initialization code...
mustart <- NULL
eval(family$initialize)
coefold <- mu <- eta <- NULL
old.pdev <- null.pdev <- sum(dev.resids(y, linkinv(null.eta), weights, theta)) + t(null.coef)%*%St%*%null.coef
# Calculating pen deviance using each initialization and choosing the best
tmp <- lapply(start, function(.st){
if (length(.st) != nvars){stop("Length of start should equal ", nvars, " and correspond to initial coefs for ", deparse(xnames))}
.eta <- offset + as.vector(if (NCOL(x) == 1) x * .st else x %*% .st)
.mu <- linkinv(.eta)
.pdev <- sum(dev.resids(y, .mu, weights, theta)) + t(.st)%*%St%*%.st
return( list("eta"=.eta, "mu"=.mu, "pdev"=.pdev, "start"=.st) )
})
tmp <- tmp[[ which.min(sapply(tmp, "[[", "pdev")) ]]
coefold <- start <- tmp$start
eta <- etaold <- tmp$eta
mu <- tmp$mu
pdev <- tmp$pdev
# BAD start (it's worse than null.coef) reset everything
if (pdev>old.pdev){
start <- coefold <- eta <- etaold <- mu <- NULL
} else { # GOOD start
old.pdev <- pdev
}
## Initialization of mu and eta (and coefold) if "start" did not do it
if(is.null(coefold)){ coefold <- null.coef }
if (is.null(eta)){
eta <- family$linkfun(mustart)
mu <- linkinv(eta)
etaold <- eta
}
conv <- boundary <- FALSE
for (iter in 1:control$maxit) { ## start of main fitting iteration
if (control$trace) cat(iter," ")
dd <- dDeta(y,mu,weights,theta,family,0) ## derivatives of deviance w.r.t. eta
# good <- is.finite(dd$Deta.Deta2)
w <- dd$Deta2 * .5;
wz <- w*(eta-offset) - .5*dd$Deta
z <- (eta-offset) - dd$Deta.Deta2
good <- is.finite(z)&is.finite(w)
if (control$trace&sum(!good)>0) cat("\n",sum(!good)," not good\n")
if (sum(!good)) {
use.wy <- TRUE
good <- is.finite(w)&is.finite(wz)
z[!is.finite(z)] <- 0 ## avoid NaN in .C call - unused anyway
} else use.wy <- family$use.wz
oo <- .C(.C_qgam_pls_fit1,
y=as.double(z[good]),X=as.double(x[good,]),w=as.double(w[good]),wy = as.double(wz[good]),
E=as.double(Sr),Es=as.double(Eb),n=as.integer(sum(good)),
q=as.integer(ncol(x)),rE=as.integer(rows.E),eta=as.double(z),
penalty=as.double(1),rank.tol=as.double(rank.tol),
nt=as.integer(control$nthreads),use.wy=as.integer(use.wy))
if (oo$n<0) { ## then problem is indefinite - switch to +ve weights for this step
if (control$trace) cat("**using positive weights\n")
# problem is that Fisher can be very poor for zeroes
## index weights that are finite and positive
good <- is.finite(dd$Deta2)
good[good] <- dd$Deta2[good]>0
w[!good] <- 0
wz <- w*(eta-offset) - .5*dd$Deta
z <- (eta-offset) - dd$Deta.Deta2
good <- is.finite(z)&is.finite(w)
if (sum(!good)) {
use.wy <- TRUE
good <- is.finite(w)&is.finite(wz)
z[!is.finite(z)] <- 0 ## avoid NaN in .C call - unused anyway
} else use.wy <- family$use.wz
oo <- .C(.C_qgam_pls_fit1, ##.C_pls_fit1,
y=as.double(z[good]),X=as.double(x[good,]),w=as.double(w[good]),wy = as.double(wz[good]),
E=as.double(Sr),Es=as.double(Eb),n=as.integer(sum(good)),
q=as.integer(ncol(x)),rE=as.integer(rows.E),eta=as.double(z),
penalty=as.double(1),rank.tol=as.double(rank.tol),
nt=as.integer(control$nthreads),use.wy=as.integer(use.wy))
}
# if(control$epsilon == Inf){ ### MATTEO ###################################
# startOld <- drop(start)
# lprOld <- drop(x%*%startOld)
# }
start <- oo$y[1:ncol(x)] ## current coefficient estimates
penalty <- oo$penalty ## size of penalty
eta <- drop(x%*%start) ## the linear predictor (less offset)
if (any(!is.finite(start))) { ## test for breakdown
conv <- FALSE
warning("Non-finite coefficients at iteration ",
iter)
return(list(REML=NA)) ## return immediately signalling failure
}
mu <- linkinv(eta <- eta + offset)
dev <- sum(dev.resids(y, mu, weights,theta))
########################## MATTEO ###################################
# if(control$epsilon == Inf)
# {
# .myObj <- function(.alpha){
# .param <- .alpha * start + (1-.alpha) * startOld
# .lpr <- .alpha * eta + (1-.alpha) * lprOld
# .mu <- linkinv(.lpr)
# .dev <- sum(dev.resids(y, .mu, weights, theta))
# return(.dev)
# }
#
# .opt <- optimize(.myObj, c(0, 2))$minimum
# cat(" ", .opt)
#
# start <- .opt * start + (1-.opt) * startOld
# eta <- .opt * eta + (1-.opt) * lprOld
# mu <- linkinv(eta)
# dev <- sum(dev.resids(y, mu, weights,theta))
#
# boundary <- TRUE
# penalty <- t(start)%*%St%*%start ## reset penalty too
# }
########################## MATTEO ###################################
## now step halve under non-finite deviance...
if (!is.finite(dev)) {
if (is.null(coefold)) {
if (is.null(null.coef))
stop("no valid set of coefficients has been found:please supply starting values",
call. = FALSE)
## Try to find feasible coefficients from the null.coef and null.eta
coefold <- null.coef
etaold <- null.eta
}
#warning("Step size truncated due to divergence",
# call. = FALSE)
ii <- 1
while (!is.finite(dev)) {
if (ii > control$maxit)
stop("inner loop 1; can't correct step size")
ii <- ii + 1
start <- (start + coefold)/2
eta <- (eta + etaold)/2
mu <- linkinv(eta)
dev <- sum(dev.resids(y, mu, weights,theta))
}
boundary <- TRUE
penalty <- t(start)%*%St%*%start ## reset penalty too
if (control$trace)
cat("Step halved: new deviance =", dev, "\n")
} ## end of infinite deviance correction
## now step halve if mu or eta are out of bounds...
if (!(valideta(eta) && validmu(mu))) {
#warning("Step size truncated: out of bounds",
# call. = FALSE)
ii <- 1
while (!(valideta(eta) && validmu(mu))) {
if (ii > control$maxit)
stop("inner loop 2; can't correct step size")
ii <- ii + 1
start <- (start + coefold)/2
eta <- (eta + etaold)/2
mu <- linkinv(eta)
}
boundary <- TRUE
dev <- sum(dev.resids(y, mu, weights))
penalty <- t(start)%*%St%*%start ## need to reset penalty too
if (control$trace)
cat("Step halved: new deviance =", dev, "\n")
} ## end of invalid mu/eta handling
## now check for divergence of penalized deviance....
pdev <- dev + penalty ## the penalized deviance
if (control$trace) cat("penalized deviance =", pdev, "\n")
div.thresh <- 10*(.1+abs(old.pdev))*.Machine$double.eps^.5
if (pdev-old.pdev>div.thresh) { ## solution diverging
ii <- 1 ## step halving counter
if (iter == 1 && (pdev-null.pdev>div.thresh)) { ## Doing worse than null.coef at 1st iterat -> shrink towards zero
etaold <- null.eta; coefold <- null.coef; old.pdev <- null.pdev
}
while (pdev - old.pdev > div.thresh) { ## step halve until pdev <= old.pdev
if (ii > 100)
stop("inner loop 3; can't correct step size")
ii <- ii + 1
start <- (start + coefold)/2
eta <- (eta + etaold)/2
mu <- linkinv(eta)
dev <- sum(dev.resids(y, mu, weights,theta))
pdev <- dev + t(start)%*%St%*%start ## the penalized deviance
if (control$trace)
cat("Step halved: new penalized deviance =", pdev, "\n")
}
} ## end of pdev divergence
## convergence testing...
if (abs(pdev - old.pdev)/(0.1 + abs(pdev)) < control$epsilon) {
## Need to check coefs converged adequately, to ensure implicit differentiation
## ok. Testing coefs unchanged is problematic under rank deficiency (not guaranteed to
## drop same parameter every iteration!)
grad <- 2 * t(x[good,])%*%((w[good]*(x%*%start)[good]-wz[good]))+ 2*St%*%start
if (max(abs(grad)) > control$epsilon*max(abs(start+coefold))/2) {
old.pdev <- pdev ## not converged quite enough
coef <- coefold <- start
etaold <- eta
##muold <- mu
} else { ## converged
conv <- TRUE
coef <- start
break
}
} else { ## not converged
old.pdev <- pdev
coef <- coefold <- start
etaold <- eta
}
} ## end of main loop
## so at this stage the model has been fully estimated
coef <- as.numeric(T %*% coef)
## now obtain derivatives, if these are needed...
# check.derivs <- FALSE
# while (check.derivs) { ## debugging code to check derivatives
# eps <- 1e-7
# fmud.test(y,mu,weights,theta,family,eps = eps)
# fetad.test(y,mu,weights,theta,family,eps = eps)
# }
##################### Calculate the Bayesian covariance matrix
if( needVb )
{
dd <- dDeta(y,mu,weights,theta,family,deriv)
w <- dd$Deta2 * .5
z <- (eta-offset) - dd$Deta.Deta2 ## - .5 * dd$Deta[good] / w
wf <- pmax(0,dd$EDeta2 * .5) ## Fisher type weights
wz <- w*(eta-offset) - 0.5*dd$Deta ## Wz finite when w==0
gdi.type <- if (any(abs(w)<.Machine$double.xmin*1e20)||any(!is.finite(z))) 1 else 0
good <- is.finite(wz)&is.finite(w)
residuals <- z - (eta - offset)
residuals[!is.finite(residuals)] <- NA
z[!is.finite(z)] <- 0 ## avoid passing NA etc to C code
ntot <- length(theta) + length(sp)
rSncol <- unlist(lapply(UrS,ncol))
## Now drop any elements of dd that have been dropped in fitting...
if (sum(!good)>0) { ## drop !good from fields of dd, weights and pseudodata
z <- z[good]; w <- w[good]; wz <- wz[good]; wf <- wf[good]
dd$Deta <- dd$Deta[good];dd$Deta2 <- dd$Deta2[good]
dd$EDeta2 <- dd$EDeta2[good]
if (deriv>0) dd$Deta3 <- dd$Deta3[good]
if (deriv>1) dd$Deta4 <- dd$Deta4[good]
if (length(theta)>1) {
if (deriv>0) {
dd$Dth <- dd$Dth[good,];
dd$Detath <- dd$Detath[good,]; dd$Deta2th <- dd$Deta2th[good,]
if (deriv>1) {
dd$Detath2 <- dd$Detath2[good,]; dd$Deta3th <- dd$Deta3th[good,]
dd$Deta2th2 <- dd$Deta2th2[good,];dd$Dth2 <- dd$Dth2[good,]
}
}
} else {
if (deriv>0) {
dd$Dth <- dd$Dth[good];
dd$Detath <- dd$Detath[good]; dd$Deta2th <- dd$Deta2th[good]
if (deriv>1) {
dd$Detath2 <- dd$Detath2[good]; dd$Deta3th <- dd$Deta3th[good]
dd$Deta2th2 <- dd$Deta2th2[good]; dd$Dth2 <- dd$Dth2[good]
}
}
}
}
oo <- .C(.C_qgam_gdi2,
X=as.double(x[good,]),E=as.double(Sr),Es=as.double(Eb),rS=as.double(unlist(rS)),
U1 = as.double(U1),sp=as.double(exp(sp)),theta=as.double(theta),
z=as.double(z),w=as.double(w),wz=as.double(wz),wf=as.double(wf),Dth=as.double(dd$Dth),
Det=as.double(dd$Deta),
Det2=as.double(dd$Deta2),Dth2=as.double(dd$Dth2),Det.th=as.double(dd$Detath),
Det2.th=as.double(dd$Deta2th),Det3=as.double(dd$Deta3),Det.th2 = as.double(dd$Detath2),
Det4 = as.double(dd$Deta4),Det3.th=as.double(dd$Deta3th), Deta2.th2=as.double(dd$Deta2th2),
beta=as.double(coef),b1=as.double(rep(0,ntot*ncol(x))),w1=as.double(rep(0,ntot*length(z))),
D1=as.double(rep(0,ntot)),D2=as.double(rep(0,ntot^2)),
P=as.double(0),P1=as.double(rep(0,ntot)),P2 = as.double(rep(0,ntot^2)),
ldet=as.double(1-2*(scoreType=="ML")),ldet1 = as.double(rep(0,ntot)),
ldet2 = as.double(rep(0,ntot^2)),
rV=as.double(rep(0,ncol(x)^2)),
rank.tol=as.double(.Machine$double.eps^.75),rank.est=as.integer(0),
n=as.integer(sum(good)),q=as.integer(ncol(x)),M=as.integer(nSp),
n.theta=as.integer(length(theta)), Mp=as.integer(Mp),Enrow=as.integer(rows.E),
rSncol=as.integer(rSncol),deriv=as.integer(deriv),
fixed.penalty = as.integer(rp$fixed.penalty),nt=as.integer(control$nthreads),
type=as.integer(gdi.type),dVkk=as.double(rep(0,nSp^2)))
rV <- matrix(oo$rV,ncol(x),ncol(x)) ## rV%*%t(rV)*scale gives covariance matrix
rV <- T %*% rV
Vb <- rV%*%t(rV)
} else {
Vb <- NULL
}
list("coefficients" = coef, "Vb" = Vb)
} ## gam.fitF
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.