Nothing
## file nlreg/R/nlreg.R, v 1.2-2.2 2019-01-30
##
## Copyright (C) 2000-2018 Ruggero Bellio & Alessandra R. Brazzale
##
## This file is part of the "nlreg" package for R. This program is
## free software; you can redistribute it and/or modify it under the
## terms of the GNU General Public License as published by the Free
## Software Foundation; either version 2 of the License, or (at your
## option) any later version.
##
## This library is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
## 02111-1307 USA or look up the web page
## http://www.gnu.org/copyleft/gpl.html.
##
## Please send any comments, suggestions or errors found to:
## Alessandra R. Brazzale, Department of Statistics, University of
## Padova, Via C. Battisti 241/243, 35121 Padova (PD), Italy.
## Email: alessandra.brazzale@unipd.it
## Web: http://www.stat.unipd.it/~brazzale
Dmean <- function(nlregObj, hessian = TRUE)
{
rc <- nlregObj$coef
do.call("deriv3", list(expr = nlregObj$meanFun, namevec = names(rc),
function.arg = c(names(nlregObj$ws$allPar), "..."),
hessian = hessian))
}
Dvar <- function(nlregObj, hessian = TRUE)
{
rc <- nlregObj$coef
vp <- nlregObj$varPar
namevec <- if(nlregObj$ws$xVar)
c(names(rc), names(vp))
else names(vp)
do.call("deriv3", list(expr = nlregObj$varFun, namevec = namevec,
function.arg = c(names(nlregObj$ws$allPar), "..."),
hessian = hessian))
}
nlreg <- function(formula, weights = NULL,
data = sys.frame(sys.parent()), start,
offset = NULL, subset = NULL,
control = list(x.tol = 1e-6, rel.tol = 1e-6,
step.min = 1/2048, maxit = 100),
trace = FALSE, hoa = FALSE)
{
m <- match.call()
if( missing(formula) )
stop("formula is missing, with no default")
if( !is.null(weights) & !is.call(weights) )
stop("variance function must be a call")
if( missing(start) )
stop("starting values are missing, with no default")
if( !is.null(offset) & is.null(names(offset)) )
stop("\'offset\' argument must be a named vector")
mFormula <- formula
mFormula[[2]] <- NULL
allx <- c(all.vars(formula), if(is.call(weights)) all.vars(weights))
if( !all(match(names(start), allx, nomatch=FALSE)) )
stop("some starting values are missing, with no default")
if( !is.null(offset) )
{
if( (names(offset) != "logs") &&
(!match(names(offset), allx, nomatch=FALSE)) )
stop("offset variable is not included in original formula")
if( match(names(offset), names(start), nomatch=FALSE) )
stop("no starting value admitted for offset parameter")
}
missingData <- missing(data)
if( !missingData )
if( !is.data.frame(data) )
stop("\"data\" is not a data frame")
if( !missing(control) )
{
if( is.null(control$x.tol) )
control$x.tol <- 1e-6
if( is.null(control$rel.tol) )
control$rel.tol <- 1e-6
if( is.null(control$step.min) )
control$step.min <- 1/2048
if( is.null(control$maxit) )
control$maxit <- 100
}
if( is.null(weights) )
{
if(trace)
cat("\nconstant variance function: switch to \"nls\"\n\n")
if( !missingData )
{
new.data <- data
if( !is.null(offset) && (names(offset) != "logs") )
{
new.data <- cbind(new.data, offset)
names(new.data)[dim(new.data)[[2]]] <- names(offset)
}
fit.temp <- list(as.name("nls"),
formula=formula, data=as.name("new.data"),
start=as.name("start"),
control=nls.control(maxiter=control$maxit,
tol=control$rel.tol,
minFactor=control$step.min),
algorithm="default", trace=trace)
if( !is.null(subset) )
fit.temp <- c(fit.temp, subset=subset)
fit.temp <- eval(as.call(fit.temp))
}
else
{
if( !is.null(offset) )
do.call("<-", list(names(offset), offset))
fit.temp <- list(as.name("nls"),
formula=formula, start=as.name("start"),
control=nls.control(maxiter=control$maxit,
tol=control$rel.tol,
minFactor=control$step.min),
algorithm="default", trace=trace)
if( !is.null(subset) )
fit.temp <- c(fit.temp, subset=subset)
fit.temp <- eval(as.call(fit.temp))
}
n <- length(residuals(fit.temp))
logs <- if( is.null(offset) || (names(offset) != "logs") )
log(sum(residuals(fit.temp)^2)/n)
else offset
if( !missingData )
{
.nd <- names(data)
.md <- match(".addVar", .nd, nomatch=0)
if(.md > 0)
{
.nd[.md] <- NA
cov.pos <- as.logical(match(.nd, all.vars(formula[[3]]), nomatch=0))
xx <- unlist(data[cov.pos])
yval <- eval(formula[[2]], data)
new.data <- list( xx=xx, repl=rep(1, n), dupl=1:n, t1=yval, t2=yval^2)
new.data <- lapply(new.data, as.vector)
names(new.data) <- c(names(data)[cov.pos], "repl", "dupl", "t1", "t2")
}
else
{
cov.pos <- match(names(data), all.vars(mFormula))
cov.pos <- !is.na(cov.pos)
new.data <- list( (xx <- unique(unlist(data[cov.pos]))),
tapply(unlist(data[cov.pos]), unlist(data[cov.pos]),
function(x) sum(duplicated(x)) + 1)[order(xx)],
match(unlist(data[cov.pos]),
unique(unlist(data[cov.pos]))),
tapply(eval(formula[[2]], data),
unlist(data[cov.pos]),
function(x) sum(x))[order(xx)],
tapply(eval(formula[[2]], data),
unlist(data[cov.pos]),
function(x) sum(x^2))[order(xx)] )
new.data <- lapply(new.data, as.vector)
names(new.data) <- c(names(data)[cov.pos], "repl", "dupl",
"t1", "t2")
}
}
else
{
yval <- eval(formula[[2]])
new.data <- list( repl=rep(1, n), dupl=1:n, t1=yval, t2=yval^2)
new.data <- lapply(new.data, as.vector)
names(new.data) <- c("repl", "dupl", "t1", "t2")
}
fit <- list( coef = coef(fit.temp),
varPar = if( is.null(offset) ||
(names(offset) != "logs") )
c(logs=logs) else NULL,
offset = offset,
logLik = -n/2*logs - n/2*log(2*pi) -
sum(residuals(fit.temp)^2)/exp(logs)/2,
meanFun = mFormula,
varFun = as.formula(call("~",
call("exp", as.name("logs")))),
data = new.data,
fitted = fitted(fit.temp),
weights = rep(exp(logs), n),
residuals = residuals(fit.temp)/sqrt(exp(logs)),
start = start,
call = m )
allpar <- if( is.null(offset) )
c(fit$coef, fit$varPar)
else c(fit$coef, fit$varPar, offset)
fit <- c(fit, list(ws = list( allPar = allpar,
homVar = TRUE,
xVar = FALSE,
hoa = hoa,
missingData = missingData,
frame = if(!missingData)
m$data else NULL ) ))
attr(fit, "class") <- c("nlreg", "nls")
if(hoa)
{
if( !any(match("nlreg.diag", unlist(lapply(sys.calls(),
function(x) x[[1]])), nomatch=FALSE)) )
cat("\ndifferentiating mean and variance function -- may take a while\n")
fit$ws$md <- Dmean(fit)
fit$ws$vd <- if( is.null(offset) || (names(offset) != "logs") )
Dvar(fit) else NULL
}
}
else
{
minRegCoef <- function(regCoef, varPar, offset, label, meanFun,
varFun, data, subset, isCoef, lastIter)
{
new.data <- data
if( is.null(new.data) )
{
new.frame <- as.list(regCoef)
names(new.frame) <- if(isCoef) label else names(regCoef)
}
else
{
class(new.data) <- "list"
tmp <- regCoef
names(tmp) <- if(isCoef) label else names(regCoef)
new.frame <- c(new.data, tmp)
}
if( !is.null(varPar) )
{
tmp <- varPar
names(tmp) <- if(isCoef) names(varPar) else label
new.frame <- c(new.frame, tmp)
}
if( !is.null(offset) )
{
tmp <- offset
names(tmp) <- names(offset)
new.frame <- c(new.frame, tmp)
}
num <- eval( call("-", meanFun[[2]], meanFun[[3]]),
envir=new.frame )
num <- num^2
den <- eval( varFun[[2]], envir=new.frame )
if( !is.null(subset) )
{
num <- num[subset] ; den <- den[subset]
}
nulloff <- is.null(offset)
if( !nulloff ) nulloff <- (names(offset) != "logs")
res <- if( nulloff )
{
if(!lastIter)
log(sum(num/den)/length(num))*length(num) +
sum(log(den)) + length(num) +
length(num)*log(2*pi)
else log(sum(num/den)/length(num))
}
else
{
sigma2 <- exp(offset)
if(!lastIter)
sum(log(den)) + sum(num/den)/sigma2 +
length(num)*log(sigma2) + length(num)*log(2*pi)
else offset
}
res
}
minVarPar <- function(varPar, regCoef, offset, label, meanFun,
varFun, data, subset)
{
res <- minRegCoef(regCoef=regCoef, varPar=varPar, offset=offset,
label=label, meanFun=meanFun, varFun=varFun,
data=data, subset=subset, isCoef=FALSE,
lastIter=FALSE)
res
}
fitted.val <- function(formula, weights, par, data, subset)
{
new.data <- data
if( is.null(data) )
{
new.frame <- as.list(par)
names(new.frame) <- names(par)
}
else
{
class(new.data) <- "list"
tmp <- as.list(par)
names(tmp) <- names(par)
new.frame <- c(new.data, tmp)
}
fitted <- eval(formula[[3]], envir=new.frame)
weights <- eval(weights[[2]], envir=new.frame)
resid <- eval( call("-", formula[[2]], formula[[3]]),
envir=new.frame )/sqrt(weights)
if( !is.null(subset) )
{
fitted <- fitted[subset]
weigths <- weights[subset]
resid <- resid[subset]
}
list(fitted = fitted, weights = weights, resid=resid)
}
para <- names(start)
coef.idx <- match(para, all.vars(formula[[3]]), nomatch=0)
var.idx <- (coef.idx == 0)
coef.idx <- (coef.idx != 0)
tol <- control$x.tol ; test1 <- 0 ; test2 <- 0
reltol <- control$rel.tol
maxit <- control$maxit
go <- TRUE
coef.old <- coef.new <- if( any(coef.idx) )
unlist(start[coef.idx]) else NULL
var.old <- var.new <- if( any(var.idx) )
unlist(start[var.idx]) else NULL
idx <- 0
if(trace) cat("\n")
while(go)
{
if( any(coef.idx) )
{
coef.new <- optim(par=coef.old, fn=minRegCoef, varPar=var.old,
offset=offset, label=names(coef.old),
meanFun=formula, varFun=weights,
subset=subset,
data=if(!missingData) data else NULL,
isCoef=TRUE, lastIter=FALSE,
# control=control, method="BFGS")
control=list(maxit=control$maxit, reltol=control$rel.tol),
method="BFGS")
log.lik <- -coef.new$value/2
coef.new <- coef.new$par
test2 <- max(abs((coef.old-coef.new)/coef.new))
coef.old <- coef.new
}
if( any(var.idx) )
{
var.new <- optim(par=var.old, fn=minVarPar, regCoef=coef.old,
offset=offset, label=names(var.old),
meanFun=formula, varFun=weights,
subset=subset,
data=if(!missingData) data else NULL,
# control=control, method="BFGS")
control=list(maxit=control$maxit, reltol=control$rel.tol),
method="BFGS")
log.lik <- -var.new$value/2
var.new <- var.new$par
test1 <- max(abs((var.old-var.new)/var.new))
var.old <- var.new
}
if( idx==0 )
{
new.lik <- log.lik
test3 <- reltol + 1
}
else
{
old.lik <- new.lik
new.lik <- log.lik
test3 <- max(abs((old.lik-new.lik)/new.lik))
}
idx <- idx + 1
if(trace)
cat(paste("iteration", idx, ": log likelihood =",
format(log.lik), "\n"))
if ( ((test1 < tol) && (test2 < tol)) || (test3 < reltol) ||
(idx == maxit) )
go <- FALSE
}
if( (maxit > 1) && (idx == maxit) )
warning(paste("\nconvergence not obtained in", maxit,
"iterations"))
coeff <- if( any(coef.idx) ) coef.old else NULL
nulloff <- is.null(offset)
if(!nulloff) nulloff <- (names(offset) != "logs")
varPar <- if( nulloff )
c(logs=minRegCoef(coef.old, var.old, offset=offset,
names(coef.old), formula, weights,
data=if(!missing(data))
data else NULL,
subset=subset, TRUE, TRUE))
else NULL
if( any(var.idx) ) varPar <- c(var.old, varPar)
allpar <- if( is.null(offset) )
c(coeff, varPar)
else c(coeff, varPar, offset)
varFun <- as.formula(call("~", call("*", weights[[2]],
call("exp", as.name("logs")))))
more.fit <- fitted.val(formula, varFun, allpar,
data=if(!missingData) data else NULL,
subset=subset)
n <- length(more.fit$resid)
if( !missing(data) )
{
.nd <- names(data)
.md <- match(".addVar", .nd, nomatch=0)
if(.md > 0)
{
.nd[.md] <- NA
cov.pos <- as.logical(match(.nd, all.vars(formula[[3]]), nomatch=0))
xx <- unlist(data[cov.pos])
yval <- eval(formula[[2]], data)
new.data <- list( xx=xx, repl=rep(1, n), 1:n, t1=yval, t2=yval^2)
new.data <- lapply(new.data, as.vector)
names(new.data) <- c(names(data)[cov.pos], "repl", "dupl", "t1", "t2")
}
else
{
cov.pos <- match(names(data), all.vars(formula[[3]])) ## 20.05.13
cov.pos <- !is.na(cov.pos)
new.data <- list( (xx <- unique(unlist(data[cov.pos]))),
tapply(unlist(data[cov.pos]), unlist(data[cov.pos]),
function(x) sum(duplicated(x)) + 1)[order(xx)],
match(unlist(data[cov.pos]),
unique(unlist(data[cov.pos]))),
tapply(eval(formula[[2]], data), unlist(data[cov.pos]),
function(x) sum(x))[order(xx)],
tapply(eval(formula[[2]], data), unlist(data[cov.pos]),
function(x) sum(x^2))[order(xx)] )
new.data <- lapply(new.data, as.vector)
names(new.data) <- c(names(data)[cov.pos], "repl", "dupl",
"t1", "t2")
}
}
else
{
yval <- eval(formula[[2]])
new.data <- list( repl=rep(1, n), 1:n, t1=yval, t2=yval^2)
new.data <- lapply(new.data, as.vector)
names(new.data) <- c("repl", "dupl", "t1", "t2")
}
xVar <- any( match(names(coeff), all.vars(weights[[2]]),
nomatch=0) )
fit <- list( coef = coeff,
varPar = varPar,
offset = offset,
logLik = - minRegCoef(coef.old, var.old,
offset=offset,
names(coef.old), formula,
weights,
data=if(!missing(data))
data
else NULL,
subset, TRUE, FALSE)/2,
meanFun = as.formula(call("~", formula[[3]])),
varFun = varFun,
data = new.data,
fitted = more.fit$fitted,
weights = more.fit$weights,
residuals = more.fit$resid,
start = start,
call = m,
ws = list( allPar = allpar,
homVar = FALSE,
xVar = xVar,
hoa = hoa,
missingData = missingData,
frame = if(!missingData)
m$data
else NULL,
iter = idx ) )
attr(fit, "class") <- c("nlreg", "nls")
if(hoa)
{
if( !any(match("nlreg.diag", unlist(lapply(sys.calls(),
function(x) x[[1]])), nomatch=FALSE)) )
cat("\ndifferentiating mean and variance function -- may take a while\n")
fit$ws$md <- Dmean(fit)
fit$ws$vd <- Dvar(fit)
}
}
fit
}
logLik.nlreg <- function(object, ...)
{
if( length(list(...)) )
warning("extra arguments discarded")
nobs <- length(object$residuals)
coef <- object$coef
varPar <- object$varPar
p <- length(coef) + length(varPar)
val <- object$logLik
attr(val, "nobs") <- nobs
attr(val, "npar") <- p
attr(val, "df") <- nobs - p
class(val) <- "logLik"
val
}
print.nlreg <- function(x, digits = max(3, getOption("digits")-3),
...)
{
if(!is.null(cl <- x$call))
{
cat("Formula:\n")
dput(x$call$formula)
cat("\nVariance function:\n")
dput(x$varFun)
}
coef <- x$coef
if( !is.null(coef) )
{
cat("\nRegression coefficients:\n")
print(coef, digits=digits, ...)
}
else cat("\nNo regression coefficient\n")
varPar <- x$varPar
if( !is.null(varPar) )
{
cat("\nVariance parameters:\n")
print(varPar, digits=digits, ...)
}
else cat("\nNo variance parameter\n")
of <- x$offset
if( !is.null(of) )
{
cat("\nParameter of interest:\n")
print(of, digits=digits, ...)
}
else cat("\nNo interest parameter\n")
nobs <- length(x$residuals)
cat("\nTotal number of observations:", nobs)
cat("\nTotal number of parameters:",
sum(length(coef)+length(varPar)))
cat("\n-2*Log Likelihood",
format(-2 * x$logLik, digits=digits), "\n")
if( x$ws$homVar )
cat("\nAlgorithm converged\n")
else cat("\nAlgorithm converged in", x$ws$iter,
if(x$ws$iter>1) "iterations\n" else "iteration\n")
invisible(x)
}
coef.nlreg <- function(object, ...)
object$coef
param <- function(object, ...)
UseMethod("param")
param.nlreg <- function(object, ...)
c(object$coef, object$varPar)
residuals.nlreg <- function(object, ...)
object$residuals
resid.nlreg <- residuals.nlreg
fitted.nlreg <- function(object, ...)
object$fitted
summary.nlreg <- function(object, observed = TRUE,
correlation = FALSE, digits = NULL, ...)
{
nlregObj <- object
rc <- nlregObj$coef
vp <- nlregObj$varPar
of <- nlregObj$offset
if( nlregObj$ws$hoa )
{
md <- nlregObj$ws$md
vd <- nlregObj$ws$vd
}
par <- nlregObj$ws$allPar
mu <- nlregObj$fitted
v <- nlregObj$weights
.probl <- ( nlregObj$ws$homVar && !is.null(of) )
if( .probl )
.probl <- ( names(of) =="logs" )
# attach(nlregObj$data, warn.conflicts = FALSE) ## 20.05.13
tmp <- as.list(par)
tmp <- c(nlregObj$data, tmp)
temp <- if( nlregObj$ws$hoa )
{
formals(md) <- tmp
do.call("md", tmp)
}
else
{
cat("\ndifferentiating mean function -- may take a while")
tmp.fun <- Dmean(nlregObj)
formals(tmp.fun) <- tmp
do.call("tmp.fun", tmp)
}
m1 <- attr(temp, "gradient") ; m1[!is.finite(m1)] <- 0
m2 <- attr(temp, "hessian") ; m2[!is.finite(m2)] <- 0
if( !.probl )
{
tmp <- as.list(par)
tmp <- c(nlregObj$data, tmp)
temp <- if( nlregObj$ws$hoa )
{
formals(vd) <- tmp
do.call("vd", tmp)
}
else
{
cat("\ndifferentiating variance function -- may take a while")
tmp.fun <- Dvar(nlregObj)
formals(tmp.fun) <- tmp
do.call("tmp.fun", tmp)
}
v1 <- attr(temp, "gradient") ; v1[!is.finite(v1)] <- 0
v2 <- attr(temp, "hessian") ; v2[!is.finite(v2)] <- 0
}
else
{
v1 <- v2 <- NULL
}
# detach(nlreg.data) ## 20.05.13
if( !nlregObj$ws$hoa ) cat("\n")
info <- if( observed )
obsInfo.nlreg(nlregObj, par, mu, v, m1, m2, v1, v2)
else expInfo.nlreg(nlregObj, par, mu, v, m1, v1)
s.e. <- sqrt(diag( (cov <- solve(qr(info))) ))
dimnames(cov) <- list(names(c(rc,vp)), names(c(rc,vp)))
idx1 <- 1:length(rc)
if( !is.null(vp) )
idx2 <- length(rc) + 1:length(vp)
coef <- matrix(0, ncol = 4, nrow = length(rc))
dimnames(coef) <- list(names(rc),
c("Estimate", "Std. Error", "z value",
"Pr(>|z|)"))
coef[, 1] <- as.vector(rc)
coef[, 2] <- s.e.[idx1]
coef[, 3] <- coef[, 1]/coef[, 2]
coef[, 4] <- 2 * pnorm( - abs(coef[, 3]) )
if( !is.null(vp) )
{
varPar <- matrix(0, nrow = length(vp), ncol=2)
dimnames(varPar) <- list(names(vp), c("Estimate", "Std. Error"))
varPar[,1] <- vp
varPar[,2] <- s.e.[idx2]
}
else varPar <- NULL
if( correlation )
{
correl <- matvec( vecmat(1/s.e., cov), 1/s.e.)
dimnames(correl) <- list(names(c(rc,vp)), names(c(rc,vp)))
}
else correl <- NULL
summary <- list( coefficients = coef,
varPar = varPar,
offset = nlregObj$offset,
residuals = nlregObj$resid,
covariance = cov,
correlation = correl,
logLik = nlregObj$logLik,
call = nlregObj$call,
digits = digits,
ws = nlregObj$ws )
attr(summary, "class") <- c("summary.nlreg", "summary.nls")
summary
}
print.summary.nlreg <- function(x,
digits = max(3, getOption("digits")-3),
signif.stars =
getOption("show.signif.stars"),
quote = TRUE, ...)
{
coef <- x$coef
varPar <- x$varPar
of <- x$offset
correl <- x$correl
if(missing(digits) && !is.null(x$digits))
digits <- x$digits
cat("Call: \n")
dput(x$call)
cat("\nRegression coefficients:\n")
printCoefmat(coef, digits=digits, signif.stars=signif.stars)
if( !is.null(varPar) )
{
cat("\nVariance parameters:\n")
print(varPar, digits=digits)
}
else cat("\nNo variance parameter\n")
if( !is.null(of) )
{
cat("\nParameter of interest:\n")
print(of, digits=digits)
}
else cat("\nNo interest parameter\n")
nobs <- length(x$residuals)
cat("\nTotal number of observations:", nobs)
npar <- length(coef[,1]) + if( is.null(varPar) )
0 else length(varPar[,1])
cat("\nTotal number of parameters:", npar)
cat("\n-2*Log Likelihood", format(-2 * x$logLik, digits=digits),
"\n")
if( x$ws$homVar )
cat("\nAlgorithm converged\n")
else cat("\nAlgorithm converged in", x$ws$iter,
if(x$ws$iter>1) "iterations\n" else "iteration\n")
if( !is.null(correl) )
{
p <- dim(correl)[2]
if(p > 1)
{
cat("\nCorrelation of parameters:\n")
ll <- lower.tri(correl)
correl[ll] <- format(round(correl[ll], digits))
correl[!ll] <- ""
print(correl[-1, - p, drop=FALSE], quote=FALSE, digits=digits)
}
}
invisible(x)
}
var2cor <- function(object, ...)
UseMethod("var2cor")
var2cor.nlreg <- function(object, ...)
{
object <- summary(object, corr=TRUE)
var2cor(object)
}
var2cor.summary.nlreg <- function(object, ...)
{
if( !is.null(object$correlation) )
object$correlation
else
{
cov <- object$covariance
s.e. <- sqrt(diag(cov))
correl <- matvec( vecmat(1/s.e., cov), 1/s.e.)
dimnames(correl) <- dimnames(object$covariance)
correl
}
}
nlreg.diag <- function(fitted, hoa = TRUE, infl = TRUE, trace = FALSE)
{
nlregObj <- fitted
rc <- nlregObj$coef
vp <- nlregObj$varPar
of <- nlregObj$offset
dupl <- nlregObj$data$dupl
par <- nlregObj$ws$allPar
mu <- nlregObj$fitted
v <- nlregObj$weights
.probl <- ( nlregObj$ws$homVar && !is.null(of) )
if(.probl)
.probl <- ( names(of) =="logs" )
followUp <- trace
wantHoa <- hoa
if(wantHoa)
{
nlregHoa <- nlregObj
nlregHoa$meanFun[[2]] <- call("+", nlregHoa$meanFun[[2]],
call("*", as.name("phi"),
as.name(".addVar")))
nlregHoa$coef <- c(nlregHoa$coef, phi=0)
nlregHoa$ws$allPar <- c(nlregHoa$ws$allPar, phi=0)
cat("\ndifferentiating mean function -- may take a while")
md <- Dmean(nlregHoa)
if(!.probl)
{
cat("\ndifferentiating variance function -- may take a while")
vd <- Dvar(nlregHoa)
}
cut.idx <- length(rc) + 1
.addVar <- rep(0, length(nlregObj$data$repl))
nlregHoa$data <- c(nlregHoa$data, list(.addVar=.addVar))
# attach(nlregHoa$data, warn.conflicts = FALSE) ## 20.05.13
new.par <- nlregHoa$ws$allPar
tmp <- as.list(new.par)
tmp <- c(nlregHoa$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.0 <- attr(temp, "gradient") ; m1.0[!is.finite(m1.0)] <- 0
m2.0 <- attr(temp, "hessian") ; m2.0[!is.finite(m2.0)] <- 0
m1.cut <- m1.0[,-cut.idx, drop=FALSE]
m2.cut <- m2.0[,-cut.idx, -cut.idx, drop=FALSE]
if(!.probl)
{
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.0 <- attr(temp, "gradient") ; v1.0[!is.finite(v1.0)] <- 0
v2.0 <- attr(temp, "hessian") ; v2.0[!is.finite(v2.0)] <- 0
if( nlregObj$ws$xVar )
{
v1.cut <- v1.0[,-cut.idx, drop=FALSE]
v2.cut <- v2.0[,-cut.idx, -cut.idx, drop=FALSE]
}
else { v1.cut <- v1.0 ; v2.cut <- v2.0 }
}
else
{
v1.0 <- v2.0 <- v1.cut <- v2.cut <- NULL
}
M1.0 <- m1.cut[dupl,]
# detach(nlregHoa$data) ## 20.05.13
i.0 <- expInfo.nlreg(nlregObj, par, mu, v, m1.cut, v1.cut)
j.0 <- obsInfo.nlreg(nlregObj, m1=m1.cut, m2=m2.cut,
v1=v1.cut, v2=v2.cut)
# if( infl && !nlregObj$ws$homVar && !is.null(nlregObj$varPar ) )
if( infl && !is.null(nlregObj$varPar ) ) ## 14.06.13
{
mle.rc <- solve(j.0[names(rc), names(rc), drop=FALSE],
j.0[names(rc), names(vp), drop=FALSE])
mle.vp <- solve(j.0[names(vp), names(vp), drop=FALSE],
j.0[names(vp), names(rc), drop=FALSE])
}
j.0 <- det(j.0)
newCall <- match.call(nlreg, nlregObj$call)
newCall$formula[[3]] <- nlregHoa$meanFun[[2]]
newStart <- as.list(newCall$start)
newStart <- c(newStart, list(phi=0))
newStart <- as.call(newStart)
newCall$start <- newStart
newCall$trace <- FALSE
if( nlregObj$ws$hoa )
newCall$hoa <- FALSE
newCall$data <- NULL
### newData <- eval(nlregObj$ws$frame) ## 24.05.13
if(!is.null(nlregObj$ws$frame))
newData <- eval(nlregObj$ws$frame)
else
{
.allVars <- all.vars(nlregObj$call$formula)
.start <- names(nlregObj$call$start)
.modVars <- .allVars[!match(.allVars, .start, nomatch=0)]
newData <- sapply(.modVars, function(x) eval(as.name(x)))
}
logLik.0 <- nlregObj$logLik
phi.h <- logLik.1 <- j.1 <- i.1 <- S.1 <- q.1 <- rep(NA, length(v))
}
else
{
# if( infl && !nlregObj$ws$homVar )
if( infl ) ## 14.06.13
{
md <- if(nlregObj$ws$hoa )
nlregObj$ws$md
else
{
cat("\ndifferentiating mean function -- may take a while")
Dmean(nlregObj)
}
vd <- if(nlregObj$ws$hoa )
nlregObj$ws$vd
else
{
cat("\ndifferentiating variance function -- may take a while")
Dvar(nlregObj)
}
}
else
{
md <- if(nlregObj$ws$hoa )
nlregObj$ws$md
else
{
cat("\ndifferentiating mean function -- may take a while")
Dmean(nlregObj, hessian=FALSE)
}
if( !.probl )
vd <- if(nlregObj$ws$hoa )
nlregObj$ws$vd
else
{
cat("\ndifferentiating variance function -- may take a while")
Dvar(nlregObj, hessian=FALSE)
}
}
# attach(nlregObj$data, warn.conflicts = FALSE) ## 20.05.13
tmp <- as.list(nlregObj$ws$allPar)
tmp <- c(nlregObj$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.0 <- attr(temp, "gradient") ; m1.0[!is.finite(m1.0)] <- 0
# if( infl && !nlregObj$ws$homVar )
if( infl ) ## 14.06.13
{
m2.0 <- attr(temp, "hessian") ; m2.0[!is.finite(m2.0)] <- 0
}
if(!.probl)
{
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.0 <- attr(temp, "gradient") ; v1.0[!is.finite(v1.0)] <- 0
# if( infl && !nlregObj$ws$homVar )
if( infl ) ## 14.06.13
{
v2.0 <- attr(temp, "hessian") ; v2.0[!is.finite(v2.0)] <- 0
}
}
else
{
v1.0 <- v2.0 <- NULL
}
M1.0 <- m1.0[dupl,]
# detach(nlregObj$data) ## 20.05.13
i.0 <- expInfo.nlreg(nlregObj, par, mu, v, m1.0, v1.0)
# if( infl && !nlregObj$ws$homVar && !is.null(nlregObj$varPar) )
if( infl && !is.null(nlregObj$varPar) ) ## 14.06.13
{
j.0 <- obsInfo.nlreg(nlregObj, m1=m1.0, m2=m2.0, v1=v1.0,
v2=v2.0)
mle.rc <- solve(j.0[names(rc), names(rc), drop=FALSE],
j.0[names(rc), names(vp), drop=FALSE])
mle.vp <- solve(j.0[names(vp), names(vp), drop=FALSE],
j.0[names(vp), names(rc), drop=FALSE])
}
}
# if( infl && !nlregObj$ws$homVar )
if( infl ) ## 14.06.13
{
newFrame <- eval(nlregObj$ws$frame)
newFrame <- c(newFrame, as.list(nlregObj$coef))
newFrame <- c(newFrame, as.list(nlregObj$varPar))
newFrame <- c(newFrame, as.list(nlregObj$offset))
logLik1 <- nlregObj$logLik
ld <- ld.rc <- ld.vp <- rep(0, length(v))
}
res <- nlregObj$resid ## standardized residuals
s.i.0 <- solve(qr(i.0))
dimnames(s.i.0) <- dimnames(i.0)
hh <- diag( M1.0 %*% s.i.0[names(rc), names(rc), drop=FALSE] %*%
t(M1.0) ) / v ## leverages
ha <- diag( M1.0 %*% solve( crossprod( M1.0, vecmat(1/v, M1.0) ),
t(M1.0) ) ) / v
## approximate leverages
r.p <- nlregObj$residuals / sqrt( 1 - hh )
## generalized Pearson residuals
r.s <- nlregObj$residuals / sqrt( 1 - ha )
## approximate studentized residuals
cook <- hh/(length(rc)*(1-hh))*r.p^2
## approximate influence for regression coefficients
if(followUp)
cat("\n")
for( sset in 1:length(v) )
{
if(wantHoa)
{
.addVar <- diag(1, length(v))[, sset]
..newData2 <- data.frame(cbind(newData, .addVar))
lv <- lm((nlregObj$residual*sqrt(nlregObj$weights))~.addVar-1)$coef
newCall$start["phi"] <- signif(lv, 4)
### attach(..newData2, warn.conflicts = FALSE) ## 24.05.13
newCall$data <- as.name("..newData2")
ctrl <- try(eval(newCall), silent=TRUE)
if( class(ctrl)[1] != "try-error" )
{
nlreg.temp <- eval(newCall)
lv <- coef(nlreg.temp)["phi"]
phi.h[sset] <- nlreg.temp$coef["phi"]
logLik.1[sset] <- nlreg.temp$logLik
# attach(nlreg.temp$data, warn.conflicts = FALSE) ## 20.05.13
tmp <- as.list(nlreg.temp$ws$allPar)
tmp <- c(nlreg.temp$data, list(.addVar=.addVar), tmp)
### tmp <- c(nlreg.temp$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.1 <- attr(temp, "gradient") ; m1.1[!is.finite(m1.1)] <- 0
m2.1 <- attr(temp, "hessian") ; m2.1[!is.finite(m2.1)] <- 0
tmp <- as.list(nlregHoa$ws$allPar)
tmp <- c(nlreg.temp$data, list(.addVar=.addVar), tmp)
### tmp <- c(nlreg.temp$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.0 <- attr(temp, "gradient") ; m1.0[!is.finite(m1.0)] <- 0
if(!.probl)
{
tmp <- as.list(nlreg.temp$ws$allPar)
tmp <- c(nlreg.temp$data, list(.addVar=.addVar), tmp)
### tmp <- c(nlreg.temp$data, tmp)
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.1 <- attr(temp, "gradient") ; v1.1[!is.finite(v1.1)] <- 0
v2.1 <- attr(temp, "hessian") ; v2.1[!is.finite(v2.1)] <- 0
tmp <- as.list(nlregHoa$ws$allPar)
tmp <- c(nlreg.temp$data, list(.addVar=.addVar), tmp)
### tmp <- c(nlreg.temp$data, tmp)
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.0 <- attr(temp, "gradient") ; v1.0[!is.finite(v1.0)] <- 0
}
else v1.1 <- v2.1 <- v1.0 <- NULL
# detach(nlreg.temp$data) ## 20.05.13
### detach(..newData2) ## 24.05.13
j.1[sset] <- det(obsInfo.nlreg(nlreg.temp, m1=m1.1, m2=m2.1,
v1=v1.1, v2=v2.1))
i.1[sset] <- det(expInfo.nlreg(nlreg.temp, m1=m1.1, v1=v1.1))
S.temp <- Shat.nlreg(nlreg.temp, nlregObj, par.0=new.par,
m1.1=m1.1, m1.0=m1.0, v1.1=v1.1, v1.0=v1.0)
q.temp <- qhat.nlreg(nlreg.temp, nlregObj, par.0=new.par,
m1.1=m1.1, v1.1=v1.1)
s.S.temp <- solve(qr(S.temp))
dimnames(s.S.temp) <- dimnames(S.temp)
q.1[sset] <- s.S.temp["phi",] %*% q.temp
S.1[sset] <- det(S.temp)
}
}
# if( infl && !nlregObj$ws$homVar ) ## 14.06.13
if( infl )
{
nlreg.temp <- update(nlregObj, subset=c(-sset), trace=FALSE)
newFrame[names(nlreg.temp$coef)] <- nlreg.temp$coef
newFrame[names(nlreg.temp$varPar)] <- nlreg.temp$varPar
num <- eval( call("-", nlregObj$call$formula[[2]],
nlregObj$call$formula[[3]]), envir=newFrame )
num <- num^2
den <- eval( nlregObj$varFun[[2]], envir=newFrame )
logLik0 <- -( sum(log(den)) + sum(num/den) +
length(num)*log(2*pi) )/2
ld[sset] <- 2/length(c(rc,vp))*(logLik1-logLik0)
if( !is.null(nlregObj$varPar) )
{
rc.cond <- rc + mle.rc %*% (vp - nlreg.temp$varPar)
vp.cond <- vp + mle.vp %*% (rc - nlreg.temp$coef)
newFrame[names(nlreg.temp$varPar)] <- vp.cond
num <- eval( call("-", nlregObj$call$formula[[2]],
nlregObj$call$formula[[3]]),
envir=newFrame )
num <- num^2
den <- eval( nlregObj$varFun[[2]], envir=newFrame )
logLik0 <- -( sum(log(den)) + sum(num/den) +
length(num)*log(2*pi) )/2
ld.rc[sset] <- 2/length(rc)*(logLik1-logLik0)
newFrame[names(nlreg.temp$coef)] <- rc.cond
newFrame[names(nlreg.temp$varPar)] <- nlreg.temp$varPar
num <- eval( call("-", nlregObj$call$formula[[2]],
nlregObj$call$formula[[3]]),
envir=newFrame )
num <- num^2
den <- eval( nlregObj$varFun[[2]], envir=newFrame )
logLik0 <- -( sum(log(den)) + sum(num/den) +
length(num)*log(2*pi) )/2
ld.vp[sset] <- 2/length(vp)*(logLik1-logLik0)
}
}
if(followUp)
{
# if(wantHoa || (infl && !nlregObj$ws$homVar))
if(wantHoa || infl) ## 14.06.13
cat(paste("\niteration", sset, "out of", length(v)))
if( wantHoa )
{
if( class(ctrl)[1] == "try-error" )
cat("\n\'nlreg\' did not converge: NAs produced in diagnotics")
}
}
else
if( wantHoa )
{
if( class(ctrl)[1] == "try-error" )
cat(paste("\n\'nlreg\' did not converge in iteration",
sset, "out of", length(v),
": NAs produced in diagnotics \n"))
}
}
if(followUp)
# if(wantHoa || (infl && !nlregObj$ws$homVar)) cat("\n")
if(wantHoa || infl) cat("\n") ## 14.06.13
if(wantHoa)
{
r.j <- sign(phi.h) * sqrt(2*(logLik.1-logLik.0))
## deletion residuals
rs.j <- r.j - log(abs(i.1*sqrt(j.0)*r.j/
sqrt(j.1)/S.1/q.1))/r.j
## r*-type residuals
}
diagObj <- list(fitted = nlregObj$fitted,
resid = res, rp = r.p, rs = r.s)
if(wantHoa)
diagObj <- c(diagObj, list( rj = r.j, rsj = rs.j ) )
diagObj <- c(diagObj, list( h = hh, ha = ha, cook = cook ) )
# if( infl && !nlregObj$ws$homVar )
if( infl ) ## 14.06.13
{
diagObj <- c(diagObj, list( ld = ld) )
if( !is.null(nlregObj$varPar) )
diagObj <- c(diagObj, list( ld.rc = ld.rc, ld.vp = ld.vp ))
}
diagObj <- c( diagObj, npar=length(rc) )
class(diagObj) <- "nlreg.diag"
diagObj
}
nlreg.diag.plots <- function(fitted, which = "all", subset = NULL,
iden = FALSE, labels = NULL, hoa = TRUE,
infl = TRUE, trace = FALSE, ret = FALSE,
...)
{
if( any(match("nlreg.diag", class(fitted), nomatch=FALSE)) )
{
plot(fitted, which=which, subset=subset, iden=iden, labels=labels,
...)
if(ret) fitted
}
else if( any(match("nlreg", class(fitted), nomatch=FALSE)) )
{
objNlreg <- nlreg.diag(fitted, hoa=hoa, infl=infl,
trace=trace)
plot(objNlreg, which=which, subset=subset, iden=iden,
labels=labels, ...)
if (ret) objNlreg
}
else
stop("missing \"nlreg\" or \"nlreg.diag\" object")
}
plot.nlreg.diag <- function(x, which = "all", subset = NULL,
iden = FALSE, labels = NULL, ...)
{
nlregdiag <- x
if( attr(nlregdiag, "class") != "nlreg.diag" )
stop("\nmissing \"nlreg.diag\" object")
if( is.null(subset) )
subset <- c( 1:length(nlregdiag$h) )
else if( is.logical(subset) )
subset <- ( 1:length(subset) )[subset]
else if( is.numeric(subset) && all(subset < 0) )
subset <- ( 1:(length(subset) +
length(nlregdiag$h)) )[subset]
else if(is.character(subset))
{
if( is.null(labels) )
labels <- subset
subset <- seq(along = subset)
}
choices <- c("Summary",
"Studentized residuals against fitted values",
"r* residuals against fitted values",
"Normal QQ-plot of studentized residuals",
"Normal QQ-plot of r* residuals",
"Cook statistic against h/(1-h)",
"Global influence against h/(1-h)",
"Cook statistic against observation number",
"Influence measures against observation number\n")
tmenu <- paste("plot:", choices)
if( which == "all" )
pick <- menu(tmenu,
title = "\n Make a plot selection (or 0 to exit)\n")
else if( !match(which, 2:length(choices), nomatch=FALSE) )
stop("choice not valid")
else pick <- which
if( pick == 0 )
stop("no graph required ! ")
old.par <- par() ; on.exit(par(old.par))
repeat
{
switch(pick,
"1" = { par(pty="s", mfrow=c(1,1))
close.screen(all.screens = TRUE)
split.screen(figs = c(2, 2))
screen(1)
## Plot the studentized residuals against the fitted values
x1 <- nlregdiag$fitted
plot(x1, nlregdiag$rp, xlab = "Fitted values",
ylab = "Studentized residuals", ...)
screen(2)
## Plot a normal QQ-plot of either the studentized residuals or the r*
## residuals
y2 <- if( !is.null(nlregdiag$rsj) ) nlregdiag$rsj
else nlregdiag$rp
x2 <- qnorm(ppoints(length(y2)))[rank(y2)]
plot(x2, y2,
xlab = paste("Quantiles of standard normal"),
ylab = paste("Ordered",
if( !is.null(nlregdiag$rsj) )
"r* residuals"
else "studentized residuals"),
xlim=c(-3,3), ylim=c(-3,3), ...)
abline(0, 1, lty = 2)
screen(3)
## Plot either the Cook statistics or the global influence measure
## against h/(1-h) and draw line to highlight possible influential and
## high leverage points.
hh <- nlregdiag$h/(1 - nlregdiag$h)
y3 <- if( is.null(nlregdiag$ld) )
nlregdiag$cook else nlregdiag$ld
plot(hh, y3, xlab = "h/(1-h)",
ylab = if( is.null(nlregdiag$ld) )
"Cook statistic"
else "Global influence", ...)
rx <- range(hh)
ry <- range(y3)
npar <- nlregdiag$npar
nobs <- length(nlregdiag$h)
cooky <- 8/(nobs - 2 * npar)
hy <- (2 * npar)/(nobs - 2 * npar)
if((cooky >= ry[1]) && (cooky <= ry[2]))
abline(h = cooky, lty = 2)
if((hy >= rx[1]) && (hy <= rx[2]))
abline(v = hy, lty = 2)
screen(4)
## Plot either the Cook statistics or the influence measures against
## the observation number in the original data set.
y4 <- if( is.null(nlregdiag$ld) ) nlregdiag$cook
else nlregdiag$ld
plot(subset, y4, xlab = "Case", type="h",
ylab = if( is.null(nlregdiag$ld) )
"Cook statistic"
else "Global influence", ...)
if((cooky >= ry[1]) && (cooky <= ry[2]))
abline(h = cooky, lty = 2)
xx <- list(x1, x2, hh, subset)
yy <- list(nlregdiag$rp, y2, y3, y4)
if(is.null(labels))
labels <- if( !is.null(names(nlregdiag$fitted)) )
names(nlregdiag$fitted)
else paste(1:length(nlregdiag$h))
yes <- iden
while(yes)
{
## If interaction with the plots is required then ask the user which
## plot they wish to interact with and then run identify() on that
## plot. When the user terminates identify(), reprompt until no
## further interaction is required and the user inputs a 0.
cat("****************************************************\n")
cat("Please Input a screen number (1,2,3 or 4)\n")
cat("0 will terminate the function \n")
num <- scan(n = 1)
check <- ((length(num) > 0) && ((num == 1) || (
num == 2) || (num == 3) || (num == 4)))
if(check)
{
cat(paste("Interactive Identification for screen",
num, "\n"))
cat("left button = Identify, center button = Exit\n")
screen(num, new = FALSE)
identify(xx[[num]], yy[[num]], labels, ...)
}
else yes <- FALSE
}
close.screen(all.screens=TRUE)
par(ask=FALSE)
},
"2" = { par(pty="s", mfrow=c(1,1))
## Plot the studentized residuals against the fitted values
x1 <- nlregdiag$fitted
plot(x1, nlregdiag$rp, xlab = "Fitted values",
ylab = "Studentized residuals", ...)
xx <- list(x1)
yy <- list(nlregdiag$rp)
},
"3" = { par(pty="s", mfrow=c(1,1))
## Plot the r* residuals (if available) against the fitted values
if( is.null(nlregdiag$rsj) )
{
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
text(0,0, "Not available")
}
else
{
x2 <- nlregdiag$fitted
plot(x2, nlregdiag$rsj, xlab = "Fitted values",
ylab = "r* residuals", ...)
xx <- list(x2)
yy <- list(nlregdiag$rsj)
}
},
"4" = { par(pty="s", mfrow=c(1,1))
## Plot a normal QQ-plot of the studentized residuals
y3 <- nlregdiag$rp
x3 <- qnorm(ppoints(length(y3)))[rank(y3)]
plot(x3, y3, xlab = "Quantiles of standard normal",
ylab = "Ordered studentized residuals",
xlim=c(-3,3), ylim=c(-3,3), ...)
abline(0, 1, lty = 2)
xx <- list(x3)
yy <- list(y3)
},
"5" = { par(pty="s", mfrow=c(1,1))
## Plot a normal QQ-plot of the r* residuals (if available)
if( is.null(nlregdiag$rsj) )
{
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
text(0,0, "Not available")
}
else
{
y4 <- nlregdiag$rsj
x4 <- qnorm(ppoints(length(y4)))[rank(y4)]
plot(x4, y4,
xlab = "Quantiles of standard normal",
ylab = "Ordered r* residuals", ylim=c(-3,3),
xlim=c(-3,3), ...)
abline(0, 1, lty = 2)
xx <- list(x4)
yy <- list(y4)
}
},
"6" = { par(pty="s", mfrow=c(1,1))
## Plot the Cook statistics against h/(1-h) and draw line to highlight
## possible influential and high leverage points.
hh <- nlregdiag$h/(1 - nlregdiag$h)
y5 <- nlregdiag$cook
plot(hh, y5, xlab = "h/(1-h)",
ylab = "Cook statistic", ...)
rx <- range(hh)
ry <- range(y5)
npar <- nlregdiag$npar
nobs <- length(nlregdiag$h)
cooky <- 8/(nobs - 2 * npar)
hy <- (2 * npar)/(nobs - 2 * npar)
if((cooky >= ry[1]) && (cooky <= ry[2]))
abline(h = cooky, lty = 2)
if((hy >= rx[1]) && (hy <= rx[2]))
abline(v = hy, lty = 2)
xx <- list(hh)
yy <- list(nlregdiag$cook)
},
"7" = { par(pty="s", mfrow=c(1,1))
## Plot the global influence (if available) against h/(1-h) and draw
## line to highlight possible influential and high leverage points.
if( is.null(nlregdiag$ld) )
{
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
text(0,0, "Not available")
}
else
{
hh <- nlregdiag$h/(1 - nlregdiag$h)
y6 <- nlregdiag$ld
plot(hh, y6, xlab = "h/(1-h)",
ylab = "Global influence", ...)
rx <- range(hh)
ry <- range(y6)
npar <- nlregdiag$npar
nobs <- length(nlregdiag$h)
cooky <- 8/(nobs - 2 * npar)
hy <- (2 * npar)/(nobs - 2 * npar)
if((cooky >= ry[1]) && (cooky <= ry[2]))
abline(h = cooky, lty = 2)
if((hy >= rx[1]) && (hy <= rx[2]))
abline(v = hy, lty = 2)
xx <- list(hh)
yy <- list(nlregdiag$ld)
}
},
"8" = { par(pty="s", mfrow=c(1,1))
## Plot the Cook statistics against the observation number in the
## original data set.
plot(subset, nlregdiag$cook, xlab = "Case",
ylab = "Cook statistic", type="h", ...)
ry <- range(nlregdiag$cook)
npar <- nlregdiag$npar
nobs <- length(nlregdiag$h)
cooky <- 8/(nobs - 2 * npar)
if((cooky >= ry[1]) && (cooky <= ry[2]))
xx <- list(subset)
yy <- list(nlregdiag$cook)
},
"9" = { par(pty="s", mfrow=c(1,1))
## Plot the influence measures (if available) against the observation
## number in the original data set.
if( is.null(nlregdiag$ld) )
{
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
text(0,0, "Not available")
}
else
{
if( is.null(nlregdiag$ld.rc) )
{
plot(subset, nlregdiag$ld, xlab = "Case",
ylab = "", main = "Global Influence",
type="h", ...)
ry <- range(nlregdiag$ld)
npar <- nlregdiag$npar
nobs <- length(nlregdiag$h)
cooky <- 8/(nobs - 2 * npar)
if((cooky >= ry[1]) && (cooky <= ry[2]))
abline(h = cooky, lty = 2)
if(is.null(labels))
labels <- if( !is.null(names(nlregdiag$fitted)) )
names(nlregdiag$fitted)
else paste(1:length(nlregdiag$h))
## If interaction with the plots is required then ask the user which
## plot they wish to interact with and then run identify() on that
## plot. When the user terminates identify(), reprompt until no
## further interaction is required and the user inputs a 0.
cat("****************************************************\n")
cat("Interactive Identification\n")
cat("left button = Identify, center button = Exit\n")
identify(subset, nlregdiag$ld, labels, ...)
}
else
{
split.screen(figs=c(1,3))
screen(1)
plot(subset, nlregdiag$ld, xlab = "Case",
ylab = "", main = "Global Influence",
type="h", ...)
ry <- range(nlregdiag$ld)
npar <- nlregdiag$npar
nobs <- length(nlregdiag$h)
cooky <- 8/(nobs - 2 * npar)
if((cooky >= ry[1]) && (cooky <= ry[2]))
abline(h = cooky, lty = 2)
screen(2)
plot(subset, nlregdiag$ld.rc, xlab = "Case",
ylab = "",
main = "Partial Influence\n(regr. coeff.)\n",
type="h", ...)
screen(3)
plot(subset, nlregdiag$ld.vp, xlab = "Case",
ylab = "",
main = "Partial Influence\n(var. param.)\n",
type="h", ...)
xx <- list(subset, subset, subset)
yy <- list(nlregdiag$ld, nlregdiag$ld.rc,
nlregdiag$ld.vp)
if(is.null(labels))
labels <- if( !is.null(names(nlregdiag$fitted)) )
names(nlregdiag$fitted)
else paste(1:length(nlregdiag$h))
yes <- iden
while(yes)
{
## If interaction with the plots is required then ask the user which
## plot they wish to interact with and then run identify() on that
## plot. When the user terminates identify(), reprompt until no
## further interaction is required and the user inputs a 0.
cat("****************************************************\n")
cat("Please Input a screen number (1,2 or 3)\n")
cat("0 will terminate the function \n")
num <- scan(n = 1)
check <- ((length(num) > 0) && ((num == 1) ||
(num == 2) || (num == 3)))
if(check)
{
cat(paste("Interactive Identification for screen",
num, "\n"))
cat("left button = Identify, center button = Exit\n")
screen(num, new = FALSE)
identify(xx[[num]], yy[[num]], labels, ...)
}
else yes <- FALSE
}
close.screen(all.screens=TRUE)
par(ask=FALSE)
}
}} )
if( (pick != 1) && (pick != 9) && iden )
{
if( !( (((pick==3) || (pick==5)) && is.null(nlregdiag$rsj)) ||
((pick==7) && is.null(nlregdiag$ld)) ) )
{
if(is.null(labels))
labels <- if( !is.null(names(nlregdiag$fitted)) )
names(nlregdiag$fitted)
else paste(1:length(nlregdiag$h))
## If interaction with the plots is required then ask the user which
## plot they wish to interact with and then run identify() on that
## plot. When the user terminates identify(), reprompt until no
## further interaction is required and the user inputs a 0.
cat("****************************************************\n")
cat("Interactive Identification\n")
cat("left button = Identify, center button = Exit\n")
identify(xx[[1]], yy[[1]], labels, ...)
}
}
if( which == "all" )
pick <- menu(tmenu,
title = "\n Make a plot selection (or 0 to exit)\n")
if( (pick == 0) || (which != "all") )
{
invisible(close.screen(all.screens=TRUE))
break
}
}
invisible()
}
profile.nlreg <- function(fitted, offset = "all", hoa = TRUE,
precision = 6, signif = 30, n = 50,
omit = 0.5, trace = FALSE, md, vd,
all = FALSE, ...)
{
m <- match.call()
nlregObj <- fitted
if( !is.null(nlregObj$offset) )
stop("\noffset parameter is already fixed!")
offsetName <- paste(substitute(offset))
if(offsetName == "all")
allProfiles.nlreg(fitted=nlregObj, hoa=hoa, signif=signif,
n=n, precision=precision, omit=omit,
trace=trace, call=m, ...)
else
{
if( !match(offsetName, names(nlregObj$ws$allPar), nomatch=FALSE) )
stop("\noffset parameter does not appear in formula")
signif <- 2*as.integer(signif/2)
followUp <- trace
wantHoa <- hoa
arg.list <- match.call(nlreg, nlregObj$call)
arg.list$trace <- FALSE
new.start <- nlregObj$ws$allPar[names(nlregObj$ws$allPar)!=
"logs"]
if( offsetName != "logs" )
{
idx <- match( offsetName, names(new.start) )
new.start <- new.start[-idx]
}
if( nlregObj$ws$homVar )
arg.list$weights <- NULL
rc <- if( (idx1 <- match(offsetName, names(nlregObj$coef),
nomatch=0)) )
nlregObj$coef[-idx1]
else nlregObj$coef
vp <- if( (idx2 <- match(offsetName, names(nlregObj$varPar),
nomatch=0)) )
nlregObj$varPar[-idx2]
else nlregObj$varPar
smaller <- c(names(rc), names(vp))
if( missing(md) )
if( nlregObj$ws$hoa )
md <- nlregObj$ws$md
else
{
cat("\ndifferentiating mean function -- may take a while")
md <- Dmean(nlregObj)
}
if( missing(vd) )
if( nlregObj$ws$hoa )
vd <- nlregObj$ws$vd
else
{
cat("\ndifferentiating variance function -- may take a while")
vd <- Dvar(nlregObj)
}
# attach(nlregObj$data, warn.conflicts = FALSE) ## 20.05.13
tmp <- as.list(nlregObj$ws$allPar)
tmp <- c(nlregObj$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.1 <- attr(temp, "gradient") ; m1.1[!is.finite(m1.1)] <- 0
m2.1 <- attr(temp, "hessian") ; m2.1[!is.finite(m2.1)] <- 0
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.1 <- attr(temp, "gradient") ; v1.1[!is.finite(v1.1)] <- 0
v2.1 <- attr(temp, "hessian") ; v2.1[!is.finite(v2.1)] <- 0
# detach(nlregObj$data) ## 20.05.13
logLik1 <- nlregObj$logLik
mle <- if( idx1 ) nlregObj$coef[idx1]
else nlregObj$varPar[idx2]
jh <- obsInfo.nlreg(nlregObj, m1=m1.1, m2=m2.1, v1=v1.1, v2=v2.1)
if(all) jh.all <- jh
s.e. <- solve(qr(jh))
dimnames(s.e.) <- dimnames(jh)
s.e. <- sqrt(s.e.[offsetName, offsetName])
if(wantHoa)
{
jh <- sqrt(det(jh))
jh.t <- jh * s.e.
ih <- expInfo.nlreg(nlregObj, m1=m1.1, v1=v1.1)
ih.t <- det(ih[smaller, smaller, drop=FALSE])
ih <- det(ih)
rh <- resid(nlregObj)
wh <- sqrt(nlregObj$weights)
db <- ncol(m1.1)
dp <- ncol(v1.1) + ifelse(!nlregObj$ws$xVar, db, 0)
npts <- nrow(m1.1)
Dvh <- matrix(0, nrow=npts, ncol=dp)
if(nlregObj$ws$xVar)
Dvh <- v1.1
else Dvh[,(db+1):dp] <- v1.1
Dmh <- m1.1
if( !nlregObj$ws$missingData )
{
Dvh <- Dvh[nlregObj$data$dupl,]
Dmh <- Dmh[nlregObj$data$dupl,]
}
Dvh <- Dvh/2/wh
Dmvh <- t(Dvh*rh)
Dmvh[1:db,] <- Dmvh[1:db,] + t(Dmh)
dimnames(Dmvh) <- list(names(nlregObj$ws$allPar), rep("", ncol(Dmvh)))
lvh <- - matrix( rowSums( t(Dvh / wh) ), ncol=1)
lthvh <- Dmvh %*% (aa <- cbind( vecmat( 1/wh^2, Dmh ) +
vecmat( rh/wh^2, Dvh[,1:db] ),
vecmat( 2*rh/wh^2, Dvh[,(db+1):dp] ) ))
# lthvh <- Dmvh %*% (aa <- cbind( vecmat( 1/wh^2, (Dmh + Dvh[,1:db]) ),
# vecmat( 2*rh/wh^2, Dvh[,(db+1):dp] ) ))
dimnames(lthvh) <- list(names(nlregObj$ws$allPar),
names(nlregObj$ws$allPar))
dimnames(aa) <- list(rep("", nrow(aa)), names(nlregObj$ws$allPar))
lchvh <- Dmvh[smaller,,drop=FALSE] %*% aa[,smaller,drop=FALSE]
if( nlregObj$ws$homVar )
lthvh[dp,1:db] <- lthvh[1:db,dp] <- 0
}
offsetVals <- c(mle-precision*s.e., mle+precision*s.e.)
offsetVals <- c(offsetVals[1],
offsetVals[1] + cumsum(rep(diff(offsetVals)/
(signif-1),(signif-1))))
offsetVals <- signif(offsetVals, 6)
x <- r <- w <- rep(NA, signif)
if(wantHoa)
{
rs.sk <- q.sk <- l.sk <- inf.sk <- np.sk <- rep(NA, signif)
rs.fr <- q.fr <- l.fr <- inf.fr <- np.fr <- rep(NA, signif)
}
mat <- matrix( rep(NA, signif*(length(nlregObj$ws$allPar) +
ifelse(wantHoa, 4, 2))),
nrow=length(offsetVals) )
new.start.0 <- new.start
go.on <- TRUE ; next.idx <- signif/2
counter <- 0 ; first.time <- TRUE ; limit <- 0
if(followUp)
if(nlregObj$ws$hoa) cat("\n")
else cat("\n\n")
for(i in 1:2)
{
while( go.on && (next.idx!=limit) )
{
of <- offsetVals[next.idx]
x[next.idx] <- mat[next.idx,1] <- of
names(of) <- offsetName
arg.list$offset <-
eval(parse(text=paste("call(\"c\",", offsetName, "=", of, ")")))
arg.list$start <- as.name("new.start.0")
if( nlregObj$ws$hoa )
arg.list$hoa <- FALSE
ctrl <- try(eval(arg.list), silent=TRUE)
if( class(ctrl)[1] != "try-error" )
{
nlregObj0 <- eval(arg.list)
logLik0 <- nlregObj0$logLik
w[next.idx] <- (mle-of)/s.e.
r[next.idx] <- sqrt( 2* (logLik1-logLik0) ) * sign(mle-of)
if(wantHoa)
{
# attach(nlregObj0$data, warn.conflicts = FALSE) ## 20.05.13
tmp <- as.list(nlregObj0$ws$allPar)
tmp <- c(nlregObj0$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.0 <- attr(temp, "gradient") ; m1.0[!is.finite(m1.0)] <- 0
m2.0 <- attr(temp, "hessian") ; m2.0[!is.finite(m2.0)] <- 0
m1.00 <- m1.0[,names(rc), drop=FALSE]
m2.00 <- m2.0[,names(rc),names(rc), drop=FALSE]
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.0 <- attr(temp, "gradient") ; v1.0[!is.finite(v1.0)] <- 0
v2.0 <- attr(temp, "hessian") ; v2.0[!is.finite(v2.0)] <- 0
if( (offsetName == "logs") && nlregObj0$ws$homVar )
{ v1.00 <- v1.0 ; v2.00 <- v2.0 }
else if( !nlregObj0$ws$xVar)
{
v1.00 <- v1.0[, names(vp), drop=FALSE]
v2.00 <- v2.0[, names(vp), names(vp), drop=FALSE]
}
else
{
v1.00 <- v1.0[, c(names(rc), names(vp)), drop=FALSE]
v2.00 <- v2.0[, c(names(rc), names(vp)),
c(names(rc), names(vp)), drop=FALSE]
}
# detach(nlregObj0$data) ## 20.05.13
jt <- obsInfo.nlreg(nlregObj0, m1=m1.00, m2=m2.00,
v1=v1.00, v2=v2.00)
Sh <- Shat.nlreg(nlregObj, nlregObj0, m1.1=m1.1, m1.0=m1.0,
v1.1=v1.1, v1.0=v1.0)
Sh.t <- Sh[smaller, smaller, drop=FALSE]
qh <- qhat.nlreg(nlregObj, nlregObj0, m1.1=m1.1, v1.1=v1.1)
s.Sh <- solve(qr(Sh))
dimnames(s.Sh) <- dimnames(Sh)
Sinvq <- s.Sh[match(offsetName,dimnames(Sh)[[1]]),] %*% qh
lambda.sk.num <- jh * det(Sh)* det(Sinvq)
lambda.sk.den <- ih * (xx <- sqrt(abs(det(jt))))
q.sk[next.idx] <- sign(r[next.idx]) *
abs(lambda.sk.num / lambda.sk.den)
lambda.sk <- abs(q.sk[next.idx] / r[next.idx])
l.sk[next.idx] <- as.numeric(lambda.sk)
np.sk[next.idx] <- log( abs(det(Sh.t) * jh.t / ih.t / xx ) ) /
r[next.idx]
inf.sk[next.idx] <- (xx <- log(l.sk[next.idx])/r[next.idx]) -
np.sk[next.idx]
rs.sk[next.idx] <- as.numeric(r[next.idx] + xx)
names(rs.sk[next.idx]) <- offsetName
rt <- resid(nlregObj0)
wt <- sqrt(nlregObj0$weights)
Dvt <- matrix(0, nrow=npts, ncol=dp)
if(nlregObj$ws$xVar)
Dvt <- v1.0
else Dvt[,(db+1):dp] <- v1.0
Dmt <- m1.0
if( !nlregObj$ws$missingData )
{
Dvt <- Dvt[nlregObj$data$dupl,]
Dmt <- Dmt[nlregObj$data$dupl,]
}
Dvt <- Dvt/2/wt
lvt <- - matrix( rowSums( matvec(Dmvh, rt/wt) ), ncol=1)
dimnames(lvt) <- list(names(nlregObj$ws$allPar), "")
if( nlregObj$ws$homVar )
lvt[names(rc),] <- 0
aa <- cbind( vecmat( 1/wt^2, Dmt ) + vecmat( rt/wt^2, Dvt[,1:db] ),
vecmat( 2*rt/wt^2, Dvt[,(db+1):dp] ) )
# aa <- cbind( vecmat( 1/wt^2, (Dmt + Dvt[,1:db]) ),
# vecmat( 2*rt/wt^2, Dvt[,(db+1):dp] ) )
lthvt <- Dmvh %*% aa
dimnames(lthvt) <- list(names(nlregObj$ws$allPar),
names(nlregObj$ws$allPar))
dimnames(aa) <- list(rep("", nrow(aa)),
names(nlregObj$ws$allPar))
lchvt <- Dmvh[smaller,,drop=FALSE] %*% aa[,smaller,drop=FALSE]
AA <- cbind(lvh-lvt, lthvt[, smaller, drop=FALSE])
lambda.fr.num <- jh * det(AA)
lambda.fr.den <- det(lthvh) * (xx <- sqrt(abs(det(jt))))
q.fr[next.idx] <- sign(r[next.idx]) *
abs(lambda.fr.num / lambda.fr.den)
lambda.fr <- abs(q.fr[next.idx] / r[next.idx])
l.fr[next.idx] <- as.numeric(lambda.fr)
np.fr[next.idx] <-
log( abs(det(lchvt[smaller, smaller, drop=FALSE]) /
det(lchvh[smaller, smaller, drop=FALSE]) *
jh.t/xx) ) / r[next.idx]
inf.fr[next.idx] <- (xx <- log(l.fr[next.idx])/r[next.idx]) -
np.fr[next.idx]
rs.fr[next.idx] <- as.numeric(r[next.idx] + xx)
names(rs.fr[next.idx]) <- offsetName
mat[next.idx,-1] <- c(nlregObj0$coef, nlregObj0$varPar,
rs.sk[next.idx], rs.fr[next.idx],
r[next.idx], w[next.idx])
go.on <- ( min(abs(c(w[next.idx], r[next.idx],
rs.sk[next.idx], rs.fr[next.idx]))) < 2.4 )
if( any(!is.finite(c(w[next.idx], r[next.idx],
rs.sk[next.idx], rs.fr[next.idx]))) )
{
x[next.idx] <- mat[next.idx,1] <- NA
w[next.idx] <- r[next.idx] <- NA
rs.sk[next.idx] <- rs.fr[next.idx] <- NA
q.sk[next.idx] <- l.sk[next.idx] <- NA
np.sk[next.idx] <- inf.sk[next.idx] <- NA
q.fr[next.idx] <- l.fr[next.idx] <- NA
np.fr[next.idx] <- inf.fr[next.idx] <- NA
go.on <- FALSE
}
else
{
next.idx <- next.idx + ifelse(first.time, -1, 1)
counter <- counter + 1
new.start.0 <- c(nlregObj0$coef,
nlregObj0$varPar[names(nlregObj0$varPar)!="logs"])
if(followUp)
{
if( !any(match("allProfiles.nlreg",
unlist(lapply(sys.calls(),
function(x) x[[1]])),
nomatch=FALSE)) )
cat(paste(m$offset, "=", format(of), "\n"))
else cat(paste(format(of), "\n"))
}
}
}
else
{
mat[next.idx,-1] <- c(nlregObj0$coef, nlregObj0$varPar,
r[next.idx], w[next.idx])
go.on <- ( min(abs(c(w[next.idx], r[next.idx]))) < 2.4 )
if( any(!is.finite(c(w[next.idx], r[next.idx]))) )
{
x[next.idx] <- mat[next.idx,1] <- NA
w[next.idx] <- r[next.idx] <- NA
go.on <- FALSE
}
else
{
next.idx <- next.idx + ifelse(first.time, -1, 1)
counter <- counter + 1
new.start.0 <- c(nlregObj0$coef,
nlregObj0$varPar[names(nlregObj0$varPar)!="logs"])
if(followUp)
{
if( !any(match("all.Profiles.nlreg",
unlist(lapply(sys.calls(),
function(x) x[[1]])),
nomatch=FALSE)) )
cat(paste(m$offset, "=", format(of), "\n"))
else cat(paste(format(of), "\n"))
}
}
}
}
else
{
go.on <- TRUE
next.idx <- next.idx + ifelse(first.time, -1, 1)
counter <- counter + 1
if(followUp)
{
if( !any(match("allProfiles.nlreg",
unlist(lapply(sys.calls(),
function(x) x[[1]])),
nomatch=FALSE)) )
cat(paste(m$offset, "=", format(of), "\n"))
else cat(paste(format(of), "\n"))
cat("\'nlreg\' did not converge: offset value not considered for profiling\n")
}
else
cat(paste("\n\'nlreg\' did not converge: offset value",
m$offset, "=", format(of),
"not considered for profiling\n"))
}
}
if(first.time)
{
if( next.idx == 0 )
cat("\nlower range too short\n")
go.on <- TRUE ; next.idx <- signif/2 + 1
new.start.0 <- new.start
first.time <- FALSE ; limit <- signif
}
else
if( next.idx == signif )
cat("\nupper range too short\n")
}
if(!followUp && nlregObj$ws$hoa) cat("\n")
x <- x[is.finite(w)] ; w <- w[is.finite(w)] ; r <- r[is.finite(r)]
if(wantHoa)
{
rs.sk <- rs.sk[is.finite(rs.sk)]
q.sk <- q.sk[is.finite(q.sk)] ; l.sk <- l.sk[is.finite(l.sk)]
np.sk <- np.sk[is.finite(np.sk)] ; inf.sk <- inf.sk[is.finite(inf.sk)]
rs.fr <- rs.fr[is.finite(rs.fr)]
q.fr <- q.fr[is.finite(q.fr)] ; l.fr <- l.fr[is.finite(l.fr)]
np.fr <- np.fr[is.finite(np.fr)] ; inf.fr <- inf.fr[is.finite(inf.fr)]
}
mat <- mat[is.finite(mat[,2]),]
dimnames(mat) <- list(c(), c(offsetName, smaller,
if(wantHoa) c("rs.sk", "rs.fr") else NULL, "r", "w"))
cond1 <- ( x < (mle-omit*s.e.) )
cond2 <- ( x > (mle+omit*s.e.) )
profile <- list(w = spline(x[cond1 | cond2], w[cond1 | cond2],
n=n),
r = spline(x[cond1 | cond2], r[cond1 | cond2],
n=n))
if(wantHoa)
profile <- c(profile,
list( rs.sk = spline(x[cond1 | cond2],
rs.sk[cond1 | cond2], n=n),
q.sk = spline(x[cond1 | cond2],
q.sk[cond1 | cond2], n=n),
np.sk = spline(x[cond1 | cond2],
np.sk[cond1 | cond2], n=n),
inf.sk = spline(x[cond1 | cond2],
inf.sk[cond1 | cond2], n=n),
rs.fr = spline(x[cond1 | cond2],
rs.fr[cond1 | cond2], n=n),
q.fr = spline(x[cond1 | cond2],
q.fr[cond1 | cond2], n=n),
np.fr = spline(x[cond1 | cond2],
np.fr[cond1 | cond2], n=n),
inf.fr = spline(x[cond1 | cond2],
inf.fr[cond1 | cond2], n=n)
))
profile <- c(profile,
list( mle = c(mle, se=s.e.),
para = mat,
obsI = if(all) jh.all else NULL,
hoa = wantHoa,
points = counter,
n = n,
call = m ))
attr(profile, "class") <- "nlreg.profile"
profile
}
}
allProfiles <- function(fitted, ...)
UseMethod("allProfiles")
allProfiles.nlreg <- function(fitted, hoa = TRUE, precision = 6,
signif = 30, n = 50, omit = 0.5,
trace = FALSE, call, ...)
{
cat("\nlong calculation --- may take a while\n")
m <- match.call()
new.call <- m
new.call[1] <- call("profile.nlreg")
nlregObj <- fitted
if( nlregObj$ws$hoa )
md <- nlregObj$ws$md
else
{
cat("\ndifferentiating mean function -- may take a while")
md <- Dmean(nlregObj)
}
if( nlregObj$ws$hoa )
vd <- nlregObj$ws$vd
else
{
cat("\ndifferentiating variance function -- may take a while")
vd <- Dvar(nlregObj)
}
new.call$md <- md
new.call$vd <- vd
new.call$all <- TRUE
followUp <- trace
wantHoa <- hoa
allOffsets <- names(nlregObj$ws$allPar)
npar <- length(allOffsets)
prof <- list()
if(followUp && !nlregObj$ws$hoa) cat("\n")
for(of in allOffsets)
{
if(followUp)
cat(paste("\nOffset:", of))
if(nlregObj$ws$hoa) cat("\n")
new.call$offset <- of
prof.temp <- eval(new.call)
prof <- c(prof,
list(list(x = prof.temp$w$x,
y = cbind(w = prof.temp$w$y,
r = prof.temp$r$y,
rs.sk = if(wantHoa) prof.temp$rs.sk$y
else NULL,
rs.fr = if(wantHoa) prof.temp$rs.fr$y
else NULL),
para = prof.temp$para,
mle = prof.temp$mle )))
}
names(prof) <- allOffsets
prof <- c(prof, list(npar=npar, obsI=prof.temp$obsI,
hoa=wantHoa, call=call))
attr(prof, "class") <- "all.nlreg.profiles"
prof
}
plot.nlreg.profile <- function(x = stop("nothing to plot"),
alpha = 0.05, add.leg = FALSE,
stats = c("sk", "fr"),
cex = 0.7, cex.lab = 1, cex.axis = 1,
cex.main = 1, lwd1 = 1, lwd2 = 2,
lty1 = "solid", lty2 = "solid",
cl1 = "blue", cl2 = "red",
col = "black", ylim = c(-3,3), ...)
{
m <- match.call()
stats <- match.arg(stats)
if( !match(stats, c("sk","fr"), nomatch=FALSE) )
stop("choice not valid for \"stats\": either \"sk\" or \"fr\"")
if(x$hoa)
cat(paste("Higher order method used:",
ifelse(stats=="sk", "Skovgaard's", "Fraser's"),
"r*\n"))
names(x)[if(stats == "sk") 3:6 else 7:10] <- c("rs", "q", "np", "inf")
confLimit <- qnorm(1-alpha/2)
plot(x$w, type="n", ylim=ylim, cex=cex, cex.lab=cex.lab,
cex.axis=cex.axis, cex.main=cex.main,
xlab = paste("parameter of interest", names(x$mle[1])),
ylab = "profiles", ...)
lines(x$w, lty="dashed", col=col, ...)
lines(x$r, lwd=lwd1, lty=lty1, col=cl1)
if(x$hoa)
lines(x$rs, lwd=lwd2, lty=lty2, col=cl2)
abline(h=confLimit, lty="dotted", col=col)
abline(h=-confLimit, lty="dotted", col=col)
if(add.leg)
{
cat("\nChoose legend position\n")
if(x$hoa)
legend(locator(1), c("Wald", "r",
ifelse(stats=="sk", "r* (Sk)", "r* (Fr)")),
lwd=c(par("lwd"), lwd1, lwd2),
lty=c("dashed", lty1, lty2), col=c(col, cl1, cl2))
else
legend(locator(1), c("Wald", "r"), lwd=c(par("lwd"), lwd1),
lty=c("dashed", lty1), col=c(col, cl1))
}
invisible(x)
}
plot.all.nlreg.profiles <- function(x = stop("nothing to plot"),
nframe, alpha = 0.05,
stats = c("sk", "fr"), cex = 0.7,
cex.lab = 1, cex.axis = 1,
cex.main = 1, lwd1 = 1, lwd2 = 2,
lty1 = "solid", lty2 = "solid",
cl1 = "blue", cl2 = "red",
col = "black", ylim = c(-3,3),
...)
{
stats <- match.arg(stats)
if( !match(stats, c("sk","fr"), nomatch=FALSE) )
stop("choice not valid for \"stats\": either \"sk\" or \"fr\"")
if(x$hoa)
cat(paste("Higher order method used:",
ifelse(stats=="sk", "Skovgaard's", "Fraser's"),
"r*\n"))
if(missing(nframe))
{
pAsk <- par("ask")
par(ask = TRUE)
on.exit( par(ask = pAsk) )
}
else
{
pMfrow <- par("mfrow")
par(mfrow = nframe)
on.exit( par(mfrow = pMfrow) )
}
npar <- x$npar
plot.prof <- function(x, alpha, stats, cex, cex.lab, cex.axis, cex.main,
lwd1, lwd2, lty1, lty2, cl1, cl2, col, ylim,
hoa, ...)
{
confLimit <- qnorm(1-alpha/2)
matplot(x$x, x$y, type="n", ylim=ylim, lty=lty1, col=col, cex=cex,
cex.lab=cex.lab, cex.axis=cex.axis, cex.main=cex.main,
xlab=paste("parameter of interest", names(x$mle[1])),
ylab="profiles", ...)
lines(x$x, x$y[,1], lty="dashed", col=col, ...)
lines(x$x, x$y[,2], lwd=lwd1, lty=lty1, col=cl1)
if(hoa)
lines(x$x, x$y[,ifelse(stats=="sk", 3, 4)], lwd=lwd2, lty=lty2,
col=cl2)
abline(h=confLimit, lty="dotted", col=col)
abline(h=-confLimit, lty="dotted", col=col)
}
sapply(x[1:npar], plot.prof, alpha=alpha, stats=stats, cex=cex,
cex.lab=cex.lab, cex.axis=cex.axis, cex.main=cex.main,
lwd1=lwd1, lwd2=lwd2, lty1=lty1, lty2=lty2, cl1=cl1, cl2=cl2,
col=col, ylim=ylim, hoa=x$hoa, ...)
invisible(x)
}
print.nlreg.profile <- function(x = stop("nothing to plot"),
alpha = 0.05, add.leg = FALSE,
stats = c("sk", "fr"),
cex = 0.7, cex.lab = 1, cex.axis = 1,
cex.main = 1, lwd1 = 1, lwd2 = 2,
lty1 = "solid", lty2 = "solid",
cl1 = "blue", cl2 = "red",
col = "black", ylim = c(-3,3), ...)
plot.nlreg.profile(x = x, alpha = alpha, add.leg = add.leg,
stats = stats, cex = cex,
cex.lab = cex.lab, cex.axis = cex.axis,
cex.main = cex.main, lwd1 = lwd1,
lwd2 = lwd2, lty1 = lty1, lty2 = lty2,
cl1 = cl1, cl2 = cl2, col = col,
ylim = ylim, ...)
print.all.nlreg.profiles <- function(x = stop("nothing to plot"),
nframe, alpha = 0.05,
stats = c("sk", "fr"), cex = 0.7,
cex.lab = 1, cex.axis = 1,
cex.main = 1, lwd1 = 1, lwd2 = 2,
lty1 = "solid", lty2 = "solid",
cl1 = "blue", cl2 = "red",
col = "black", ylim = c(-3,3),
...)
plot.all.nlreg.profiles(x = x, nframe = nframe, alpha = alpha,
stats = stats, cex = cex,
cex.lab = cex.lab, cex.axis = cex.axis,
cex.main = cex.main, lwd1 = lwd1,
lwd2 = lwd2, lty1 = lty1, lty2 = lty2,
cl1 = cl1, cl2 = cl2, col = col,
ylim = ylim, ...)
contour.all.nlreg.profiles <-
function(x, offset1, offset2, alpha = c(0.1, 0.05),
stats = c("sk", "fr"), ret = FALSE,
plotit = TRUE, drawlabels = FALSE, lwd1 = 1, lwd2 = 1,
lty1 = "solid", lty2 = "solid", cl1 = "blue",
cl2 = "red", col = "black", pch1 = 1, pch2 = 16,
cex = 0.5, ...)
{
m <- match.call()
plotIt <- plotit
if( !plotIt && !ret )
stop("There is nothing to be returned ...")
stats <- match.arg(stats)
if( !match(stats, c("sk", "fr"), nomatch=FALSE) )
stop("choice not valid for \"stats\": either \"sk\" or \"fr\"")
which <- ifelse( stats=="sk", "rs.sk", "rs.fr" )
profilesObj <- x
hoa <- profilesObj$hoa
if(hoa)
cat(paste("Higher order method used:",
ifelse(stats=="sk", "Skovgaard's", "Fraser's"),
"r*\n"))
npar <- profilesObj$npar
notGiven <- ( missing(offset1) && missing(offset2) )
if( !notGiven )
{
if( (missing(offset1) && !missing(offset2)) ||
(!missing(offset1) && missing(offset2)) )
stop("specify either no or both interest parameters")
para1 <- paste(substitute(offset1))
para2 <- paste(substitute(offset2))
if( para1 == para2 )
stop("interest parameters must be different")
para <- c(para1, para2)
.switch <- match(para, names(profilesObj)[1:npar])
if( .switch[2] < .switch[1] )
para <- c(para2, para1)
npar <- 2
}
else para <- names(profilesObj)[1:npar]
if(plotIt)
{
pPty <- par("pty")
par(pty = "s")
on.exit( par(pty = pPty) )
if( notGiven )
{
split.screen(figs = c(npar,npar))
on.exit( close.screen(all.screens = TRUE) )
}
else
{
pMfrow <- par("mfrow")
par(mfrow = c(1,2))
on.exit( par(mfrow = pMfrow) )
}
}
if(ret)
contourObj <- list()
mle.cov <- solve(qr(profilesObj$obsI))
dimnames(mle.cov) <- dimnames(profilesObj$obsI)
x.w0.r <- seq(-3, 3, length=10)
x.w.r <- as.matrix(expand.grid(x.w0.r, x.w0.r))
quantChi <- sqrt(qchisq(1-alpha, 2))
quantNorm <- qnorm(1-alpha/2)
pos1 <- pos2 <- idx <- 0
for( para1 in para)
{
pos1 <- pos1 + 1
obj1 <- profilesObj[[para1]]
mle1 <- obj1$mle
x.w1 <- seq(-3*mle1[2], 3*mle1[2], length=10)
if( notGiven && plotIt )
{
screen(pos1+npar*pos2)
plot(obj1$x, obj1$y[,"w"], type="n", xlab="", ylab="",
main=para1, ylim=c(-3,3), ...)
lines(obj1$x, obj1$y[,"w"], lty="dashed", col=col)
lines(obj1$x, obj1$y[,"r"], lwd=lwd1, lty=lty1, col=cl1)
if(hoa)
lines(obj1$x, obj1$y[,which], lwd=lwd2, lty=lty2, col=cl2)
abline(h=quantNorm, lty="dotted", col=col)
abline(h=-quantNorm, lty="dotted", col=col)
}
if( ret && notGiven )
contourObj <- c(contourObj,
list(list(profile=list(x=obj1$x, y=obj1$y))))
obj1 <- obj1$para
if( pos1 != npar )
for( para2 in para[(pos1+1):npar] )
{
idx <- idx + 1
obj2 <- profilesObj[[para2]]
mle2 <- obj2$mle
x.w2 <- seq(-3*mle2[2], 3*mle2[2], length=10)
x.w <- as.matrix(expand.grid(x.w1, x.w2))
a.cov <- mle.cov[c(para1,para2), c(para1,para2)]
a.corr <- sqrt(diag(a.cov))
a.corr <- vecmat( 1/a.corr, matvec(a.cov, 1/a.corr) )
val.w <- matrix(diag(x.w %*% solve(a.cov, t(x.w))),
ncol=10)
val.wr <- matrix(diag(x.w.r %*% solve(a.corr, t(x.w.r))),
ncol=10)
obj2 <- obj2$para
para2.r1 <- smooth.spline(obj1[,"r"], obj1[,para2],
all.knots=TRUE)
r2.para2 <- smooth.spline(obj2[,para2], obj2[,"r"],
all.knots=TRUE)
para2.rs1 <- smooth.spline(obj1[,which], obj1[,para2],
all.knots=TRUE)
rs2.para2 <- smooth.spline(obj2[,para2], obj2[,which],
all.knots=TRUE)
tau.r1 <- tau.r2 <- para.r1 <- para.r2 <- list()
tau.rs1 <- tau.rs2 <- para.rs1 <- para.rs2 <- list()
avg.r <- avg.rs <- list()
for(qq in seq(along=quantChi))
{
cut.r1 <- c(-1, 1)
temp <- predict(para2.r1, c(-1,1)*quantChi[qq])$y
temp <- predict(r2.para2, temp)$y
cut.r2 <- c(temp/quantChi[qq], -1, 1)
para1.r2 <- smooth.spline(obj2[,"r"], obj2[,para1],
all.knots=TRUE)
r1.para1 <- smooth.spline(obj1[,para1], obj1[,"r"],
all.knots=TRUE)
temp <- predict(para1.r2, c(-1,1)*quantChi[qq])$y
temp <- predict(r1.para1, temp)$y
cut.r1 <- c(cut.r1, temp/quantChi[qq])
para1.r1 <- smooth.spline(obj1[,"r"], obj1[,para1],
all.knots=TRUE)
para2.r2 <- smooth.spline(obj2[,"r"], obj2[,para2],
all.knots=TRUE)
avg.r <- c(avg.r, list((acos(cut.r1) + acos(cut.r2))/2))
diff.r <- acos(cut.r1) - acos(cut.r2)
avg.r[[qq]][(diff.r<0) & !is.na(diff.r)] <-
-avg.r[[qq]][(diff.r<0) & !is.na(diff.r)]
diff.r <- abs(diff.r)
if( !any(is.na(avg.r[[qq]])) )
{
inter.r <- smooth.spline(avg.r[[qq]], diff.r,
all.knots=TRUE)
new.diff.r <- predict(inter.r, pi)$y
avg.r[[qq]] <- c(-pi, avg.r[[qq]], pi)
diff.r <- c(new.diff.r, diff.r, new.diff.r)
inter.r <- spline(avg.r[[qq]], diff.r, 300)
tau.r1 <- c(tau.r1, list(quantChi[qq] *
cos(inter.r$x+inter.r$y/2)))
tau.r2 <- c(tau.r2, list(quantChi[qq] *
cos(inter.r$x-inter.r$y/2)))
para.r1 <- c(para.r1, list(predict(para1.r1,
tau.r1[[qq]])$y))
para.r2 <- c(para.r2, list(predict(para2.r2,
tau.r2[[qq]])$y))
}
else
{
cut1.tmp <- cut.r1*quantChi[qq]
cut2.tmp <- cut.r2*quantChi[qq]
tau.r1 <- c(tau.r1, list(cut1.tmp))
tau.r2 <- c(tau.r2, list(cut2.tmp))
temp.1 <- predict(para1.r1, cut1.tmp)$y
temp.2 <- predict(para2.r2, cut2.tmp)$y
para.r1 <- c(para.r1, list(temp.1))
para.r2 <- c(para.r2, list(temp.2))
}
if(hoa)
{
cut.rs1 <- c(-1, 1)
temp <- predict(para2.rs1, c(-1,1)*quantChi[qq])$y
temp <- predict(rs2.para2, temp)$y
cut.rs2 <- c(temp/quantChi[qq], -1, 1)
para1.rs2 <- smooth.spline(obj2[,which], obj2[,para1],
all.knots=TRUE)
rs1.para1 <- smooth.spline(obj1[,para1], obj1[,which],
all.knots=TRUE)
temp <- predict(para1.rs2, c(-1,1)*quantChi[qq])$y
temp <- predict(rs1.para1, temp)$y
cut.rs1 <- c(cut.rs1, temp/quantChi[qq])
avg.rs <- c(avg.rs, list((acos(cut.rs1) +
acos(cut.rs2))/2))
diff.rs <- acos(cut.rs1) - acos(cut.rs2)
avg.rs[[qq]][(diff.rs<0) & !is.na(diff.rs)] <-
-avg.rs[[qq]][(diff.rs<0) & !is.na(diff.rs)]
diff.rs <- abs(diff.rs)
if( !any(is.na(avg.rs[[qq]])) )
{
inter.rs <- smooth.spline(avg.rs[[qq]], diff.rs,
all.knots=TRUE)
new.diff.rs <- predict(inter.rs, pi)$y
avg.rs[[qq]] <- c(-pi, avg.rs[[qq]], pi)
diff.rs <- c(new.diff.rs, diff.rs, new.diff.rs)
inter.rs <- spline(avg.rs[[qq]], diff.rs, 300)
tau.rs1 <- c(tau.rs1, list(quantChi[qq] *
cos(inter.rs$x+inter.rs$y/2)))
tau.rs2 <- c(tau.rs2, list(quantChi[qq] *
cos(inter.rs$x-inter.rs$y/2)))
r1.rs1 <- smooth.spline(obj1[,which], obj1[,"r"],
all.knots=TRUE)
r2.rs2 <- smooth.spline(obj2[,which], obj2[,"r"],
all.knots=TRUE)
tau.rs1[[qq]] <- predict(r1.rs1, tau.rs1[[qq]])$y
tau.rs2[[qq]] <- predict(r2.rs2, tau.rs2[[qq]])$y
para.rs1 <- c(para.rs1, list(predict(para1.r1,
tau.rs1[[qq]])$y))
para.rs2 <- c(para.rs2, list(predict(para2.r2,
tau.rs2[[qq]])$y))
}
else
{
cut1.tmp <- cut.rs1*quantChi[qq]
cut2.tmp <- cut.rs2*quantChi[qq]
r1.rs1 <- smooth.spline(obj1[,which], obj1[,"r"],
all.knots=TRUE)
r2.rs2 <- smooth.spline(obj2[,which], obj2[,"r"],
all.knots=TRUE)
temp.1 <- predict(r1.rs1, cut1.tmp)$y
temp.2 <- predict(r2.rs2, cut2.tmp)$y
tau.rs1 <- c(tau.rs1, list(temp.1))
tau.rs2 <- c(tau.rs2, list(temp.2))
temp.1 <- predict(para1.r1, temp.1)$y
temp.2 <- predict(para2.r2, temp.2)$y
para.rs1 <- c(para.rs1, list(temp.1))
para.rs2 <- c(para.rs2, list(temp.2))
}
}
}
if( notGiven && plotIt )
screen(pos1+idx+npar*pos2)
xlim <- c(min(c(unlist(tau.r1),
if(hoa) unlist(tau.rs1) else NULL)) ,
max(c(unlist(tau.r1),
if(hoa) unlist(tau.rs1) else NULL)))
xlim <- xlim + c(-0.1,0.1)*diff(xlim)
ylim <- c(min(c(unlist(tau.r2),
if(hoa) unlist(tau.rs2) else NULL)) ,
max(c(unlist(tau.r2),
if(hoa) unlist(tau.rs2) else NULL)))
ylim <- ylim + c(-0.1,0.1)*diff(ylim)
if( plotIt )
{
if( notGiven )
plot(tau.r1[[1]], tau.r2[[1]], xlim=xlim, ylim=ylim,
xlab = "", ylab = "", type="n", ...)
else
plot(tau.r1[[1]], tau.r2[[1]], xlim=xlim, ylim=ylim,
xlab = para1, ylab = para2, type="n",
main="N(0,1) scale", ...)
for(qq in seq(along=quantChi))
{
if( !any(is.na(avg.r[[qq]])) )
lines(tau.r1[[qq]], tau.r2[[qq]], lwd=lwd1, lty=lty1,
col=cl1)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para1, "and", para2, "!"))
points(tau.r1[[qq]], tau.r2[[qq]], pch=pch1, col=cl1,
cex=cex)
}
if( hoa )
{
if( !any(is.na(avg.rs[[qq]])) )
lines(tau.rs1[[qq]], tau.rs2[[qq]], lwd=lwd2,
lty=lty2, col=cl2)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para1, "and", para2, "!"))
points(tau.rs1[[qq]], tau.rs2[[qq]], pch=pch2,
col=cl2, cex=cex)
}
}
}
contour(x.w0.r, x.w0.r, val.wr, levels=qchisq(1-alpha, 2),
add=TRUE, lty="dashed", col=col,
labels=paste(1-alpha), drawlabels=drawlabels, ...)
}
if(ret)
w.r <- list(x.w0.r=x.w0.r, x.w0.r=x.w0.r, val.wr=val.wr,
levels=qchisq(1-alpha,2))
temp1 <- predict(para2.r1, seq(xlim[1], xlim[2], length=50))
temp2 <- predict(r2.para2, temp1$y)
trace1.r <- list(x=temp1$x, y=temp2$y)
temp1 <- predict(para1.r2, seq(ylim[1], ylim[2], length=50))
temp2 <- predict(r1.para1, temp1$y)
trace2.r <- list(x=temp2$y, y=temp1$x)
if( plotIt )
{
lines(trace1.r, lty="solid", col=col)
lines(trace2.r, lty="solid", col=col)
}
if( notGiven && plotIt )
screen(pos1+npar*(pos2+idx))
xlim <- c(min(c(unlist(para.r1),
if(hoa) unlist(para.rs1) else NULL)) ,
max(c(unlist(para.r1),
if(hoa) unlist(para.rs1) else NULL)))
xlim <- xlim + c(-0.1,0.1)*diff(xlim)
ylim <- c(min(c(unlist(para.r2),
if(hoa) unlist(para.rs2) else NULL)) ,
max(c(unlist(para.r2),
if(hoa) unlist(para.rs2) else NULL)))
ylim <- ylim + c(-0.1,0.1)*diff(ylim)
if( plotIt )
{
if( notGiven )
plot(para.r1[[1]], para.r2[[1]], xlim=xlim, ylim=ylim,
xlab = "", ylab = "", type="n", ...)
else
plot(para.r1[[1]], para.r2[[1]], xlim=xlim, ylim=ylim,
xlab = para1, ylab = para2, type="n",
main = "original scale", ...)
for(qq in seq(along=quantChi))
{
if( !any(is.na(avg.r[[qq]])) )
lines(para.r1[[qq]], para.r2[[qq]], lwd=lwd1, lty=lty1,
col=cl1)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para1, "and", para2, "!"))
points(para.r1[[qq]], para.r2[[qq]], pch=pch1, col=cl1,
cex=cex)
}
if(hoa)
{
if( !any(is.na(avg.rs[[qq]])) )
lines(para.rs1[[qq]], para.rs2[[qq]], lwd=lwd2,
lty=lty2, col=cl2)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para1, "and", para2, "!"))
points(para.rs1[[qq]], para.rs2[[qq]], pch=pch2,
col=cl2, cex=cex)
}
}
}
contour(x.w1+mle1[1], x.w2+mle2[1], val.w,
levels=qchisq(1-alpha, 2), add=TRUE, lty="dashed",
col=col, labels=paste(1-alpha),
drawlabels=drawlabels, ...)
}
if(ret)
w.para <- list(x.w1.mle=x.w1+mle1[1], x.w2.mle=x.w2+mle2[1],
val.w=val.w, levels=qchisq(1-alpha, 2))
trace1 <- predict(smooth.spline(obj1[,para1], obj1[,para2],
all.knots=TRUE),
seq(xlim[1], xlim[2], length=50))
temp <- predict(smooth.spline(obj2[,para2], obj2[,para1],
all.knots=TRUE), seq(ylim[1], ylim[2],
length=50))
trace2 <- list(x=temp$y, y=temp$x)
if( plotIt )
{
lines(trace1, lty="solid", col=col)
lines(trace2, lty="solid", col=col)
}
if(ret)
if( notGiven )
contourObj[[length(contourObj)]] <-
c(contourObj[[length(contourObj)]],
list(list( r = list(x=tau.r1, y=tau.r2),
rs = if(hoa)
list(x=tau.rs1, y=tau.rs2)
else NULL,
w = w.r,
trace1.r = trace1.r,
trace2.r = trace2.r,
para.r = list(x=para.r1, y=para.r2),
para.rs = if(hoa)
list(x=para.rs1, y=para.rs2)
else NULL,
para.w = w.para,
trace1.par = trace1,
trace2.par = trace2)))
else
contourObj <- c(contourObj,
list( r = list(x=tau.r1, y=tau.r2),
rs = if(hoa)
list(x=tau.rs1, y=tau.rs2)
else NULL,
w = w.r,
trace1.r = trace1.r,
trace2.r = trace2.r,
para.r = list(x=para.r1, y=para.r2),
para.rs = if(hoa)
list(x=para.rs1, y=para.rs2)
else NULL,
para.w = w.para,
trace1.par = trace1,
trace2.par = trace2,
para = c(para1, para2) ))
}
if(ret && notGiven)
names(contourObj[[length(contourObj)]]) <-
c("profile", if(pos1!=npar) para[(pos1+1):npar]
else NULL)
pos2 <- pos2 + 1
idx <- 0
}
if(ret)
{
if( notGiven )
names(contourObj) <- para
contourObj <- c(contourObj, list(npar=npar, alpha=alpha, hoa=hoa,
stats=stats, fixed=!notGiven, call=m))
attr(contourObj, "class") <- "nlreg.contours"
contourObj
}
else invisible(profilesObj)
}
plot.nlreg.contours <- function(x, alpha = c(0.1, 0.05),
drawlabels = FALSE, lwd1 = 1,
lwd2 = 1, lty1 = "solid",
lty2 = "solid", cl1 = "blue",
cl2 = "red", col = "black", pch1 = 1,
pch2 = 16, cex = 0.5, ...)
{
hoa <- x$hoa
npar <- x$npar
fixed <- x$fixed
alpha <- x$alpha
stats <- x$stats
which <- ifelse( stats=="sk", "rs.sk", "rs.fr" )
quantChi <- sqrt(qchisq(1-alpha, 2))
quantNorm <- qnorm(1-alpha/2)
if(hoa)
cat(paste("Higher order method used:",
ifelse(stats=="sk", "Skovgaard's", "Fraser's"),
"r*\n"))
pPty <- par("pty")
par(pty="s")
on.exit( par(pty = pPty) )
if( !fixed )
{
split.screen(figs = c(npar,npar))
on.exit( close.screen(all.screens = TRUE) )
para <- names(x)[1:npar]
pos1 <- pos2 <- idx <- 0
for( para1 in para )
{
pos1 <- pos1 + 1
obj1 <- x[[para1]]
screen(pos1+npar*pos2)
plot(obj1$profile$x, obj1$profile$y[,"w"], type="n", xlab = "",
ylab = "", main = para1, ylim=c(-3,3), ...)
lines(obj1$profile$x, obj1$profile$y[,"w"], lty="dashed",
col=col)
lines(obj1$profile$x, obj1$profile$y[,"r"], lwd=lwd1, lty=lty1,
col=cl1)
if(hoa)
lines(obj1$profile$x, obj1$profile$y[,which], lwd=lwd2,
lty=lty2, col=cl2)
abline(h=quantNorm, lty="dotted", col=col)
abline(h=-quantNorm, lty="dotted", col=col)
if( pos1 != npar )
for( para2 in para[(pos1+1):npar] )
{
idx <- idx + 1
obj2 <- obj1[[para2]]
screen(pos1+idx+npar*pos2)
xlim <- c(min(c(unlist(obj2$r$x), unlist(obj2$rs$x))) ,
max(c(unlist(obj2$r$x), unlist(obj2$rs$x))))
xlim <- xlim + c(-0.1,0.1)*diff(xlim)
ylim <- c(min(c(unlist(obj2$r$y), unlist(obj2$rs$y))) ,
max(c(unlist(obj2$r$y), unlist(obj2$rs$y))))
ylim <- ylim + c(-0.1,0.1)*diff(ylim)
plot(obj2$trace1.r, xlab = "", ylab = "", type="l",
xlim=xlim, ylim=ylim, col=col, ...)
lines(obj2$trace2.r, lty="solid", col=col)
for(qq in 1:length(obj2$r$x))
{
if( length(obj2$r$x[[qq]]) > 4 )
lines(obj2$r$x[[qq]], obj2$r$y[[qq]], lwd=lwd1,
lty=lty1, col=cl1)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para1, "and", para2, "!"))
points(obj2$r$x[[qq]], obj2$r$y[[qq]], pch=pch1,
col=cl1, cex=cex)
}
if( hoa )
if( length(obj2$rs$x[[qq]]) > 4 )
lines(obj2$rs$x[[qq]], obj2$rs$y[[qq]], lwd=lwd2,
lty=lty2, col=cl2)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para1, "and", para2, "!"))
points(obj2$rs$x[[qq]], obj2$rs$y[[qq]], pch=pch2,
col=cl2, cex=cex)
}
}
contour(obj2$w$x.w0.r, obj2$w$x.w0.r, obj2$w$val.wr,
levels=obj2$w$levels, lty="dashed", col=col,
add=TRUE, labels=paste(1-alpha),
drawlabels=drawlabels, ...)
screen(pos1+npar*(pos2+idx))
xlim <- c(min(c(unlist(obj2$para.r$x),
unlist(obj2$para.rs$x))) ,
max(c(unlist(obj2$para.r$x),
unlist(obj2$para.rs$x))))
xlim <- xlim + c(-0.1,0.1)*diff(xlim)
ylim <- c(min(c(unlist(obj2$para.r$y),
unlist(obj2$para.rs$y))) ,
max(c(unlist(obj2$para.r$y),
unlist(obj2$para.rs$y))))
ylim <- ylim + c(-0.1,0.1)*diff(ylim)
plot(obj2$trace1.par$x, obj2$trace1.par$y, xlab = "",
ylab = "", xlim=xlim, ylim=ylim, type="l", col=col,
...)
lines(obj2$trace2.par$x, obj2$trace2.par$y, lty="solid",
col=col)
for(qq in 1:length(obj2$para.r$x))
{
if( length(obj2$para.r$x[[qq]]) > 4 )
lines(obj2$para.r$x[[qq]], obj2$para.r$y[[qq]],
lwd=lwd1, col=cl1, lty=lty1)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para1, "and", para2, "!"))
points(obj2$para.r$x[[qq]], obj2$para.r$y[[qq]],
pch=pch1, col=cl1, cex=cex)
}
if( hoa )
if( length(obj2$para.rs$x[[qq]]) > 4 )
lines(obj2$para.rs$x[[qq]], obj2$para.rs$y[[qq]],
lwd=lwd2, lty=lty2, col=cl2)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para1, "and", para2, "!"))
points(obj2$para.rs$x[[qq]], obj2$para.rs$y[[qq]],
pch=pch2, col=cl2, cex=cex)
}
}
contour(obj2$para.w$x.w1.mle, obj2$para.w$x.w2.mle,
obj2$para.w$val.w, levels=obj2$para.w$levels,
lty="dashed", col=col, add=TRUE,
labels=paste(1-alpha), drawlabels=drawlabels,
...)
}
pos2 <- pos2 + 1
idx <- 0
}
}
else
{
para <- x$para
pMfrow <- par("mfrow")
par(mfrow=c(1,2))
on.exit( par(mfrow = pMfrow) )
xlim <- c(min(c(unlist(x$r$x), unlist(x$rs$x))) ,
max(c(unlist(x$r$x), unlist(x$rs$x))))
xlim <- xlim + c(-0.1,0.1)*diff(xlim)
ylim <- c(min(c(unlist(x$r$y), unlist(x$rs$y))) ,
max(c(unlist(x$r$y), unlist(x$rs$y))))
ylim <- ylim + c(-0.1,0.1)*diff(ylim)
plot(x$trace1.r, xlab = x$para[1], ylab = x$para[2],
main = "N(0,1) scale", type="l", xlim=xlim, ylim=ylim,
col=col, ...)
lines(x$trace2.r, lty="solid", col=col, ...)
for(qq in 1:length(x$r$x))
{
if( length(x$r$x[[qq]]) > 4 )
lines(x$r$x[[qq]], x$r$y[[qq]], lwd=lwd1, lty=lty1, col=cl1)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para[1], "and", para[2], "!"))
points(x$r$x[[qq]], x$r$y[[qq]], pch=pch1, col=cl1, cex=cex)
}
if( hoa )
if( length(x$rs$x[[qq]]) > 4 )
lines(x$rs$x[[qq]], x$rs$y[[qq]], lwd=lwd2, lty=lty2,
col=cl2)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para[1], "and", para[2], "!"))
points(x$rs$x[[qq]], x$rs$y[[qq]], pch=pch2, col=cl2,
cex=cex)
}
}
contour(x$w$x.w0.r, x$w$x.w0.r, x$w$val.wr, levels=x$w$levels,
lty="dashed", col=col, add=TRUE,
labels=paste(1-alpha), drawlabels=drawlabels, ...)
xlim <- c(min(c(unlist(x$para.r$x), unlist(x$para.rs$x))) ,
max(c(unlist(x$para.r$x), unlist(x$para.rs$x))))
xlim <- xlim + c(-0.1,0.1)*diff(xlim)
ylim <- c(min(c(unlist(x$para.r$y), unlist(x$para.rs$y))) ,
max(c(unlist(x$para.r$y), unlist(x$para.rs$y))))
ylim <- ylim + c(-0.1,0.1)*diff(ylim)
plot(x$trace1.par, xlab = x$para[1], ylab = x$para[2],
main = "original scale", type="l", xlim=xlim, ylim=ylim,
col=col, ...)
lines(x$trace2.par, lty="solid", col=col)
for(qq in 1:length(x$para.r$x))
{
if( length(x$para.r$x[[qq]]) > 4 )
lines(x$para.r$x[[qq]], x$para.r$y[[qq]], lwd=lwd1, lty=lty1,
col=cl1)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para[1], "and", para[2], "!"))
points(x$para.r$x[[qq]], x$para.r$y[[qq]], pch=pch1, col=cl1,
cex=cex)
}
if( hoa )
if( length(x$para.rs$x[[qq]]) > 4 )
lines(x$para.rs$x[[qq]], x$para.rs$y[[qq]], lwd=lwd2,
lty=lty2, col=cl2)
else
{
warning(paste("could not calculate profile pairs sketches for parameters",
para[1], "and", para[2], "!"))
points(x$para.rs$x[[qq]], x$para.rs$y[[qq]], pch=pch2,
col=cl2, cex=cex)
}
}
contour(x$para.w$x.w1.mle, x$para.w$x.w2.mle, x$para.w$val.w,
levels=x$para.w$levels, lty="dashed", col=col,
add=TRUE, labels=paste(1-alpha), drawlabels=drawlabels,
...)
}
invisible(x)
}
print.nlreg.contours <- function(x, alpha = c(0.1, 0.05),
drawlabels = FALSE, lwd1 = 1,
lwd2 = 1, lty1 = "solid",
lty2 = "solid", cl1 = "blue",
cl2 = "red", col = "black", pch1 = 1,
pch2 = 16, cex = 0.5, ...)
plot.nlreg.contours(x = x, alpha = alpha,
drawlabels = drawlabels, lwd1 = lwd1,
lwd2 = lwd2, lty1 = lty1, lty2 = lty2,
cl1 = cl1, cl2 = cl2, col = col,
pch1 = pch1, pch2 = pch2, cex = cex, ...)
summary.nlreg.profile <- function(object, alpha = 0.05, twoside = TRUE,
digits = NULL, ...)
{
m <- match.call()
profile.obj <- object
quant <- if( twoside ) 1-alpha/2
else 1-alpha
quant <- qnorm(1-quant) ; quant <- c(quant, -quant)
wald.ci <- matrix(profile.obj$mle[1] + profile.obj$mle[2]*quant,
ncol=2)
c1 <- ifelse(wald.ci[,1] < wald.ci[,2], wald.ci[,1], wald.ci[,2])
c2 <- ifelse(wald.ci[,1] < wald.ci[,2], wald.ci[,2], wald.ci[,1])
wald.ci <- cbind(c1, c2)
dimnames(wald.ci) <- list(paste("Wald (", 1-alpha, ")", sep=""),
c("lower", "upper"))
r.ci <- predict(smooth.spline(profile.obj$r$y, profile.obj$r$x),
quant)
r.ci <- matrix(r.ci$y, ncol=2)
c1 <- ifelse(r.ci[,1] < r.ci[,2], r.ci[,1], r.ci[,2])
c2 <- ifelse(r.ci[,1] < r.ci[,2], r.ci[,2], r.ci[,1])
r.ci <- cbind(c1, c2)
dimnames(r.ci) <- list(paste("r (", 1-alpha, ")", sep=""),
c("lower", "upper"))
if(profile.obj$hoa)
{
rs.sk.ci <- predict(smooth.spline(profile.obj$rs.sk$y,
profile.obj$rs.sk$x), quant)
rs.sk.ci <- matrix(rs.sk.ci$y, ncol=2)
c1 <- ifelse(rs.sk.ci[,1] < rs.sk.ci[,2], rs.sk.ci[,1], rs.sk.ci[,2])
c2 <- ifelse(rs.sk.ci[,1] < rs.sk.ci[,2], rs.sk.ci[,2], rs.sk.ci[,1])
rs.sk.ci <- cbind(c1, c2)
dimnames(rs.sk.ci) <- list(paste("r* - Sk (", 1-alpha, ")", sep=""),
c("lower", "upper"))
rs.fr.ci <- predict(smooth.spline(profile.obj$rs.fr$y,
profile.obj$rs.fr$x), quant)
rs.fr.ci <- matrix(rs.fr.ci$y, ncol=2)
c1 <- ifelse(rs.fr.ci[,1] < rs.fr.ci[,2], rs.fr.ci[,1], rs.fr.ci[,2])
c2 <- ifelse(rs.fr.ci[,1] < rs.fr.ci[,2], rs.fr.ci[,2], rs.fr.ci[,1])
rs.fr.ci <- cbind(c1, c2)
dimnames(rs.fr.ci) <- list(paste("r* - Fr (", 1-alpha, ")", sep=""),
c("lower", "upper"))
mat <- rbind(rs.fr.ci, rs.sk.ci, r.ci, wald.ci)
}
else mat <- rbind(r.ci, wald.ci)
summary.obj <- list( CI = mat )
if(profile.obj$hoa)
summary.obj <- c( summary.obj,
list( inf.sk = max(abs(profile.obj$inf.sk$y)),
np.sk = max(abs(profile.obj$np.sk$y)),
inf.fr = max(abs(profile.obj$inf.fr$y)),
np.fr = max(abs(profile.obj$np.fr$y)) ) )
mle <- profile.obj$mle
offset <- names(mle)[1]
mle <- matrix(mle, nrow=1)
dimnames(mle) <- list(offset, c("Estimate", "Std. Error"))
summary.obj <- c(summary.obj,
list(mle = mle, offset = offset, alpha = alpha,
twoside = twoside, points = profile.obj$points,
n = profile.obj$n, hoa = profile.obj$hoa,
digits = digits, call = m ))
attr(summary.obj, "class") <- "summary.nlreg.profile"
summary.obj
}
summary.all.nlreg.profiles <- function(object, alpha = 0.05,
twoside = TRUE, digits = NULL,
...)
{
m <- match.call()
profile.obj <- object
quant <- if( twoside ) 1-alpha/2
else 1-alpha
quant <- qnorm(1-quant) ; quant <- c(quant, -quant)
npar <- profile.obj$npar
getCI <- function(x, quant, alpha, hoa, ...)
{
wald.ci <- matrix(x$mle[1] + x$mle[2]*quant, ncol=2)
c1 <- ifelse(wald.ci[,1] < wald.ci[,2], wald.ci[,1], wald.ci[,2])
c2 <- ifelse(wald.ci[,1] < wald.ci[,2], wald.ci[,2], wald.ci[,1])
wald.ci <- cbind(c1, c2)
dimnames(wald.ci) <- list(paste("Wald (", 1-alpha, ")", sep=""),
c("lower", "upper"))
r.ci <- predict(smooth.spline(x$y[,2], x$x), quant)
r.ci <- matrix(r.ci$y, ncol=2)
c1 <- ifelse(r.ci[,1] < r.ci[,2], r.ci[,1], r.ci[,2])
c2 <- ifelse(r.ci[,1] < r.ci[,2], r.ci[,2], r.ci[,1])
r.ci <- cbind(c1, c2)
dimnames(r.ci) <- list(paste("r (", 1-alpha, ")", sep=""),
c("lower", "upper"))
if(hoa)
{
rs.sk.ci <- predict(smooth.spline(x$y[,3], x$x), quant)
rs.sk.ci <- matrix(rs.sk.ci$y, ncol=2)
c1 <- ifelse(rs.sk.ci[,1] < rs.sk.ci[,2], rs.sk.ci[,1], rs.sk.ci[,2])
c2 <- ifelse(rs.sk.ci[,1] < rs.sk.ci[,2], rs.sk.ci[,2], rs.sk.ci[,1])
rs.sk.ci <- cbind(c1, c2)
dimnames(rs.sk.ci) <- list(paste("r* - Sk (", 1-alpha, ")", sep=""),
c("lower", "upper"))
rs.fr.ci <- predict(smooth.spline(x$y[,4], x$x), quant)
rs.fr.ci <- matrix(rs.fr.ci$y, ncol=2)
c1 <- ifelse(rs.fr.ci[,1] < rs.fr.ci[,2], rs.fr.ci[,1], rs.fr.ci[,2])
c2 <- ifelse(rs.fr.ci[,1] < rs.fr.ci[,2], rs.fr.ci[,2], rs.fr.ci[,1])
rs.fr.ci <- cbind(c1, c2)
dimnames(rs.fr.ci) <- list(paste("r* - Fr (", 1-alpha, ")", sep=""),
c("lower", "upper"))
mat <- rbind(rs.fr.ci, rs.sk.ci, r.ci, wald.ci)
}
else mat <- rbind(r.ci, wald.ci)
mat
}
summary.obj <- lapply(profile.obj[1:npar], getCI, alpha=alpha,
quant=quant, hoa=profile.obj$hoa)
mle <- matrix(unlist(lapply(profile.obj[1:npar], "[[", "mle")),
nrow=2)
offset <- names(profile.obj)[1:npar]
dimnames(mle) <- list(c("Estimate", "Std. Error"), offset)
summary.obj <- c(summary.obj,
list(mle= mle, offset = offset, alpha = alpha,
twoside = twoside,
hoa = profile.obj$hoa, digits = digits,
call = m ))
class(summary.obj) <- "summary.all.nlreg.profiles"
summary.obj
}
print.summary.nlreg.profile <- function(x,
digits = if(!is.null(x$digits)) x$digits
else max(3, getOption("digits")-3), ...)
{
cat(paste("\n", if(x$twoside) "Two-sided" else "One-sided",
"confidence intervals for", x$offset, "\n"))
nalpha <- length(x$alpha)
idx <- rep(nalpha*(0:4), nalpha) + rep(1:nalpha, rep(5,nalpha))
xCI <- rbind(x$CI, matrix(rep(NA, 2*nalpha), ncol=2))
print(xCI[idx,], digits=digits, na.print="")
print(x$mle, digits=digits)
cat("\n")
cat(paste(x$points, "points calculated exactly\n"))
cat(paste(x$n, "points used in spline interpolation\n"))
if(x$hoa)
{
cat(paste("\nINF (Sk):", signif(x$inf.sk, digits=digits)))
cat(paste("\nINF (Fr):", signif(x$inf.fr, digits=digits)))
cat(paste("\nNP (Sk):", signif(x$np.sk, digits=digits)))
cat(paste("\nNP (Fr):", signif(x$np.fr, digits=digits), "\n"))
}
}
print.summary.all.nlreg.profiles <-
function(x, digits = if(!is.null(x$digits)) x$digits
else max(3, getOption("digits")-3), ...)
{
par.names <- x$offset
nalpha <- length(x$alpha)
idxa <- rep(nalpha*(0:4), nalpha) + rep(1:nalpha, rep(5,nalpha))
cat(paste("\n", if(x$twoside) "Two-sided" else "One-sided",
"confidence intervals for \n\n"))
for(idx in 1:length(par.names))
{
cat(paste(par.names[idx], ":",
signif(x$mle[1,idx], digits=digits), "(",
signif(x$mle[2,idx], digits=digits), ")", "\n"))
xCI <- rbind(x[[idx]], matrix(rep(NA, 2*nalpha), ncol=2))
print(xCI[idxa,], digits=digits, na.print="")
}
invisible(x)
}
mpl <- function(fitted, ...)
UseMethod("mpl")
mpl.nlreg <- function(fitted, offset = NULL, stats = c("sk", "fr"),
control = list(x.tol = 1e-6, rel.tol=1e-6,
step.min=1/2048, maxit = 100),
trace = FALSE, ...)
{
m <- match.call()
nlregObj <- fitted
if( missing(nlregObj) )
stop("fitted model object is missing, with no default")
if( nlregObj$ws$homVar )
stop("models with constant variance are not considered")
if( !is.null(nlregObj$offset) )
stop("offset parameter cannot be fixed")
if( !is.null(offset) )
if( length(offset) != length(nlregObj$varPar) )
stop("variance parameters must either be fixed or free")
stats <- match.arg(stats)
if( !match(stats, c("sk","fr"), nomatch=FALSE) )
stop("choice not valid for \"stats\": either \"sk\" or \"fr\"")
if( !is.null(offset) )
offset <- offset[match(names(offset), names(nlregObj$varPar))]
if( !missing(control) )
{
if( is.null(control$x.tol) )
control$x.tol <- 1e-6
if( is.null(control$rel.tol) )
control$rel.tol <- 1e-6
if( is.null(control$step.min) )
control$step.min <- 1/2048
if( is.null(control$maxit) )
control$maxit <- 100
}
n <- length(nlregObj$fitted)
md <- if( nlregObj$ws$hoa ) nlregObj$ws$md
else
{
cat("\ndifferentiating mean function -- may take a while")
Dmean(nlregObj)
}
vd <- if( nlregObj$ws$hoa ) nlregObj$ws$vd
else
{
cat("\ndifferentiating variance function -- may take a while")
Dvar(nlregObj)
}
# attach(nlregObj$data, warn.conflicts = FALSE) ## 20.05.13
tmp <- as.list(nlregObj$ws$allPar)
tmp <- c(nlregObj$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1 <- attr(temp, "gradient") ; m1[!is.finite(m1)] <- 0
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1 <- attr(temp, "gradient") ; v1[!is.finite(v1)] <- 0
# detach(nlregObj$data) ## 20.05.13
if( stats == "fr" )
{
rh <- resid(nlregObj)
wh <- sqrt(nlregObj$weights)
db <- ncol(m1)
dp <- ncol(v1) + ifelse(!nlregObj$ws$xVar, db, 0)
npts <- nrow(m1)
Dvh <- matrix(0, nrow=npts, ncol=dp)
if(nlregObj$ws$xVar)
Dvh <- v1
else Dvh[,(db+1):dp] <- v1
Dmh <- m1
if( !nlregObj$ws$missingData )
{
Dvh <- Dvh[nlregObj$data$dupl,]
Dmh <- Dmh[nlregObj$data$dupl,]
}
Dvh <- Dvh/2/wh
Dmvh <- t(Dvh*rh)
Dmvh[1:db,] <- Dmvh[1:db,] + t(Dmh)
}
minRegCoef <- function(regCoef, varPar, md, vd, m1, v1, nlregObj,
stats, lastIter=FALSE, is.varPar=FALSE)
{
nlreg.temp <- nlregObj
new.data <- eval(nlreg.temp$call$data)
if( is.null(new.data) )
{
newFrame <- as.list(regCoef)
names(newFrame) <- names(nlreg.temp$coef)
}
else
{
class(new.data) <- "list"
tmp <- as.list(regCoef)
names(tmp) <- names(nlreg.temp$coef)
newFrame <- c(new.data, tmp)
}
tmp <- as.list(varPar)
names(tmp) <- names(nlreg.temp$varPar)
newFrame <- c(newFrame, tmp)
meanFun <- nlreg.temp$call$formula
varFun <- nlreg.temp$varFun
num <- eval( call("-", meanFun[[2]], meanFun[[3]]),
envir=newFrame )
num <- num^2
den <- eval( varFun[[2]], envir=newFrame )
logLik <- sum(log(den)) + sum(num/den)
if( is.varPar )
{
rc <- names(nlregObj$coef) ; vp <- names(nlregObj$varPar)
if( stats == "fr" )
dimnames(Dmvh) <- list(c(rc, vp), rep("", ncol(Dmvh)))
nlreg.temp$coef <- regCoef ; names(nlreg.temp$coef) <- rc
nlreg.temp$varPar <- varPar ; names(nlreg.temp$varPar) <- vp
nlreg.temp$fitted <- eval( meanFun[[3]], envir=newFrame )
nlreg.temp$weights <- den
nlreg.temp$ws$allPar <- c(regCoef, varPar)
names(nlreg.temp$ws$allPar) <- c(rc, vp)
# attach(nlreg.temp$data, warn.conflicts = FALSE) ## 20.05.13
tmp <- as.list(nlreg.temp$ws$allPar)
tmp <- c(nlreg.temp$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.t <- attr(temp, "gradient") ; m1.t[!is.finite(m1.t)] <- 0
m2.t <- attr(temp, "hessian") ; m2.t[!is.finite(m2.t)] <- 0
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.t <- attr(temp, "gradient") ; v1.t[!is.finite(v1.t)] <- 0
v2.t <- attr(temp, "hessian") ; v2.t[!is.finite(v2.t)] <- 0
# detach(nlreg.temp$data) ## 20.05.13
if( stats == "sk" )
{
S.hat <- Shat.nlreg(nlregObj, nlreg.temp, m1.1=m1, m1.0=m1.t,
v1.1=v1, v1.0=v1.t)
S.hat <- S.hat[rc, rc]
}
else
{
rt <- sqrt(num/den)
wt <- sqrt(den)
Dvt <- matrix(0, nrow=npts, ncol=dp)
if(nlregObj$ws$xVar)
Dvt <- v1.t
else Dvt[,(db+1):dp] <- v1.t
Dmt <- m1.t
if( !nlregObj$ws$missingData )
{
Dvt <- Dvt[nlregObj$data$dupl,]
Dmt <- Dmt[nlregObj$data$dupl,]
}
Dvt <- Dvt/2/wt
lvt <- - matrix( rowSums( matvec(Dmvh, rt/wt) ), ncol=1)
dimnames(lvt) <- list(names(nlregObj$ws$allPar), "")
aa <- cbind( vecmat( 1/wt^2, Dmt ) + vecmat( rt/wt^2, Dvt[,1:db] ),
vecmat( 2*rt/wt^2, Dvt[,(db+1):dp] ) )
# aa <- cbind( vecmat( 1/wt^2, (Dmt + Dvt[,1:db]) ),
# vecmat( 2*rt/wt^2, Dvt[,(db+1):dp] ) )
lchvt <- Dmvh[rc,,drop=FALSE] %*% aa[,rc,drop=FALSE]
}
i.obs <- obsInfo.nlreg(nlreg.temp, m1=m1.t, m2=m2.t,
v1=v1.t, v2=v2.t)
if( lastIter )
{
s.i.obs <- solve(qr(i.obs))
dimnames(s.i.obs) <- dimnames(i.obs)
cov.var <- s.i.obs[vp, vp]
}
i.obs <- i.obs[rc, rc]
if( lastIter )
{
cov.coef <- solve(qr(i.obs))
dimnames(cov.coef) <- list(rc, rc)
}
lmp <- logLik - log(det(i.obs)) +
2*log(det( if(stats == "sk") S.hat else lchvt ))
# lmp <- logLik - determinant(i.obs)$mod +
# 2*determinant( if(stats == "sk") S.hat else lchvt )$mod
if( !lastIter )
lmp
else
{
mod.tmp <- list(cov.var = cov.var, cov.coef = cov.coef,
lmp = -as.vector(lmp)/2 - n/2*log(2*pi),
lp = -logLik/2 - n/2*log(2*pi) )
mod.tmp
}
}
else logLik
}
minVarPar <- function(varPar, regCoef, md, vd, m1, v1, nlregObj, stats)
{
lmp <- minRegCoef(regCoef=regCoef, varPar=varPar, md=md, vd=vd,
m1=m1, v1=v1, nlregObj=nlregObj, stats=stats,
is.varPar=TRUE)
lmp
}
tol <- control$x.tol ; test1 <- 0 ; test2 <- 0 ; test3 <- 0
reltol <- control$rel.tol
maxit <- control$maxit
go <- TRUE
coef.old <- coef.new <- nlregObj$coef
var.old <- var.new <- if( is.null(offset) ) nlregObj$varPar
else offset
idx <- 0
if(trace)
if(nlregObj$ws$hoa) cat("\n")
else cat("\n\n")
while(go)
{
if( is.null(offset) )
{
var.new <- optim(par=var.old, fn=minVarPar, regCoef=coef.old,
md=md, vd=vd, m1=m1, v1=v1, nlregObj=nlregObj,
stats=stats,
# control=control, method="BFGS")
control=list(maxit=control$maxit, reltol=control$rel.tol),
method="BFGS")
var.new <- var.new$p
test1 <- max(abs((var.old-var.new)/var.new))
var.old <- var.new
}
coef.new <- optim(par=coef.old, fn=minRegCoef, varPar=var.new,
md=md, vd=vd, m1=m1, v1=v1, nlregObj=nlregObj,
stats=stats,
# control=control, method="BFGS")
control=list(maxit=control$maxit, reltol=control$rel.tol),
method="BFGS")
if( !is.null(offset) )
lmp <- -coef.new$value/2 - n/2*log(2*pi)
coef.new <- coef.new$p
test2 <- max(abs((coef.old-coef.new)/coef.new))
coef.old <- coef.new
if( is.null(offset) )
{
lmp <- minRegCoef(coef.new, var.new, md, vd, m1, v1, nlregObj,
stats=stats, lastIter=TRUE, is.varPar=TRUE)$lmp
if( idx==0 )
{
lmp.new <- lmp
test3 <- reltol + 1
}
else
{
lmp.old <- lmp.new
lmp.new <- lmp
test3 <- max(abs((lmp.old-lmp.new)/lmp.new))
}
}
if( is.null(offset) ) idx <- idx + 1
if( trace && is.null(offset) )
cat(paste("iteration", idx,
": modified profile log likelihood =", format(lmp),
"\n"))
if ( ((test1 < tol) && (test2 < tol)) || (test3 < reltol) ||
(idx == maxit) )
go <- FALSE
}
if( (maxit>1) && (idx==maxit) )
warning(paste("\nlinear convergence not obtained in", maxit,
"iterations"))
mod.temp <- minRegCoef(coef.new, var.new, md, vd, m1, v1, nlregObj,
stats=stats, lastIter=TRUE, is.varPar=TRUE)
mplObj <- list( varPar = if( is.null(offset) ) var.new
else NULL,
coefficients = coef.new,
offset = offset,
varParMLE = nlregObj$varPar,
coefMLE = nlregObj$coef,
varParCov = if( is.null(offset) )
mod.temp$cov.var else NULL,
coefCov = mod.temp$cov.coef,
lmp = mod.temp$lmp,
lp = mod.temp$lp,
stats = stats,
formula = nlregObj$call$formula,
meanFun = nlregObj$meanFun,
varFun = nlregObj$varFun,
data = nlregObj$data,
nobs = length(nlregObj$residuals),
iter = if(is.null(offset)) idx else NULL,
call = m,
ws = nlregObj$ws )
mplObj$ws$iter <- idx
attr(mplObj, "class") <- c("mpl", "nlreg")
mplObj
}
print.mpl <- function(x, digits = max(3, getOption("digits")-3), ...)
{
stats <- x$stats
if( !is.null(cl <- x$call) )
{
cat("Formula:\n")
dput(x$formula)
cat("Variance function:\n")
dput(x$varFun)
}
cat(paste("\nHigher order method used:",
ifelse(stats=="sk", "Skovgaard's", "Fraser's"),
"r*\n"))
varPar <- x$varPar
vp <- names(varPar)
if( !is.null(varPar) )
{
varPar <- cbind(varPar, x$varParMLE)
dimnames(varPar) <- list(vp, c("MMPLE ", "MLE "))
cat("\nVariance parameters\n")
print(varPar, digits=digits, ...)
}
else
{
cat("\nVariance parameters fixed to:\n")
print(x$offset, digits=digits, ...)
}
coef <- x$coefficients
rc <- names(coef)
if( !is.null(coef) )
{
coef <- cbind(coef, x$coefMLE)
dimnames(coef) <- list(rc, c("MMPLE ", "MLE "))
cat("\nRegression coefficients\n")
print(coef, digits=digits, ...)
}
else cat("\nNo regression coefficient\n")
cat("\nTotal number of observations:", x$nobs)
cat("\nTotal number of parameters:", sum(length(rc)+length(vp)))
cat("\n-2*Log Lmp", format(-2 * x$lmp, digits=digits), "\n")
if( !is.null(x$iter) )
cat("\nAlgorithm converged in", x$iter,
if(x$iter > 1) "iterations\n" else "iteration\n")
invisible(x)
}
summary.mpl <- function(object, correlation = FALSE, digits = NULL,
...)
{
mplObj <- object
mplObj$ws$corr <- correlation
mplObj$digits <- digits
class(mplObj) <- "summary.mpl"
mplObj
}
var2cor.mpl <- function(object, ...)
{
object <- summary(object)
var2cor(object)
}
var2cor.summary.mpl <- function(object, ...)
{
if( !is.null(object$varPar) )
{
s.e. <- sqrt(diag( cov <- object$varParCov ))
corrVarPar <- vecmat( 1/s.e., matvec(cov, 1/s.e.) )
dimnames(corrVarPar) <- dimnames(object$varParCov)
}
else corrVarPar <- NULL
if( !is.null(object$coefficients) )
{
s.e. <- sqrt(diag( cov <- object$coefCov ))
corrCoef <- vecmat( 1/s.e., matvec(cov, 1/s.e.) )
dimnames(corrCoef) <- dimnames(object$coefCov)
}
else corrCoef <- NULL
list( corrVarPar=corrVarPar, corrCoef=corrCoef )
}
print.summary.mpl <- function(x, corr = FALSE,
digits = if(!is.null(x$digits)) x$digits
else max(3, getOption("digits")-3),
...)
{
stats <- x$stats
if(!is.null(cl <- x$call))
{
cat("Formula:\n")
dput(x$formula)
cat("Variance function:\n")
dput(x$varFun)
}
cat(paste("\nHigher order method used:",
ifelse(stats=="sk", "Skovgaard's", "Fraser's"),
"r*\n"))
varPar <- x$varPar
vp <- names(varPar)
if( !is.null(varPar) )
{
var.se <- sqrt(diag(x$varParCov))
varPar <- cbind(varPar, x$varParMLE, var.se)
dimnames(varPar) <- list(vp, c("MMPLE", "MLE", "Std. Error"))
cat("\nVariance parameters\n")
print(varPar, digits=digits, ...)
}
else
{
cat("\nVariance parameters fixed to:\n")
print(x$offset, digits=digits, ...)
}
coef <- x$coefficients
rc <- names(coef)
if( !is.null(coef) )
{
coef.se <- sqrt(diag(x$coefCov))
coef <- cbind(coef, x$coefMLE, coef.se)
dimnames(coef) <- list(rc, c("MMPLE", "MLE", "Std. Error"))
cat("\nRegression coefficients\n")
print(coef, digits=digits, ...)
}
else cat("\nNo regression coefficient\n")
cat("\nTotal number of observations:", x$nobs)
cat("\nTotal number of parameters:", sum(length(rc)+length(vp)))
cat("\n-2*Log Lmp", format(-2 * x$lmp, digits=digits), "\n")
if( missing(corr) )
corr <- x$ws$corr
if( corr )
{
if( !is.null(varPar) )
{
cat("\nCorrelation of variance parameters\n")
p <- length(vp)
corr <- vecmat( 1/var.se, matvec(x$varParCov, 1/var.se) )
dimnames(corr) <- dimnames(x$varParCov)
ll <- lower.tri(corr)
corr[ll] <- format(round(corr[ll], digits))
corr[!ll] <- ""
print(corr[-1, -p, drop=FALSE], quote=FALSE, digits=digits)
}
if( !is.null(rc) )
{
p <- length(rc)
if(p > 1)
{
cat("\nCorrelation of regression coefficients\n")
corr <- vecmat( 1/coef.se, matvec(x$coefCov, 1/coef.se) )
dimnames(corr) <- dimnames(x$coefCov)
ll <- lower.tri(corr)
corr[ll] <- format(round(corr[ll], digits))
corr[!ll] <- ""
print(corr[-1, -p, drop=FALSE], quote=FALSE, digits=digits)
}
}
}
invisible(x)
}
obsInfo <- function(object, ...)
UseMethod("obsInfo")
obsInfo.nlreg <- function(object, par, mu, v, m1 = NULL, m2 = NULL,
v1 = NULL, v2 = NULL, ...)
{
nlregObj <- object
rc <- nlregObj$coef
vp <- nlregObj$varPar
of <- nlregObj$offset
t1 <- nlregObj$data$t1
t2 <- nlregObj$data$t2
repl <- nlregObj$data$repl
if( nlregObj$ws$hoa )
{
md <- nlregObj$ws$md
vd <- nlregObj$ws$vd
}
if( missing(par) ) par <- nlregObj$ws$allPar
if( missing(mu) ) mu <- nlregObj$fitted
if( missing(v) ) v <- nlregObj$weights
.probl <- ( nlregObj$ws$homVar && !is.null(of) )
if( .probl )
.probl <- ( names(of) =="logs" )
# attach(nlregObj$data, warn.conflicts = FALSE) ## 20.05.13
# on.exit( detach(nlregObj$data) )
if( !nlregObj$ws$missingData )
{
mu <- mu[!duplicated(nlregObj$data$dupl)]
v <- v[!duplicated(nlregObj$data$dupl)]
}
if( is.null(m1) || is.null(m2) )
{
tmp <- as.list(par)
tmp <- c(nlregObj$data, tmp)
temp <- if( nlregObj$ws$hoa )
{
formals(md) <- tmp
do.call("md", tmp)
}
else
{
tmp.fun <- Dmean(nlregObj)
formals(tmp.fun) <- tmp
do.call("tmp.fun", tmp)
}
m1 <- attr(temp, "gradient") ; m1[!is.finite(m1)] <- 0
m2 <- attr(temp, "hessian") ; m2[!is.finite(m2)] <- 0
}
if( !.probl && ( is.null(v1) || is.null(v2) ) )
{
tmp <- as.list(par)
tmp <- c(nlregObj$data, tmp)
temp <- if( nlregObj$ws$hoa )
{
formals(vd) <- tmp
do.call("vd", tmp)
}
else
{
tmp.fun <- Dvar(nlregObj)
formals(tmp.fun) <- tmp
do.call("tmp.fun", tmp)
}
v1 <- attr(temp, "gradient") ; v1[!is.finite(v1)] <- 0
v2 <- attr(temp, "hessian") ; v2[!is.finite(v2)] <- 0
}
if( nlregObj$ws$homVar )
{
v <- v[1]
if( !.probl )
{
v1 <- v1[1,1]
v2 <- v2[1,1,1]
}
}
si <- t2 + mu^2*repl - 2*t1*mu
sir <- 2*m1*(mu*repl - t1)
if( .probl )
{
mat <- array(0, dim=rep(length(rc), 2))
mat <- colSums( ((mu*repl-t1)/v)*m2[,,,drop=FALSE] ) +
crossprod( m1, vecmat(repl/v, m1) )
}
else
{
mat <- array(0, dim=rep(length(c(rc,vp)), 2))
idx1 <- 1:length(nlregObj$coef)
if( !is.null(vp) ) idx2 <- length(idx1) + 1:length(vp)
mat[idx1,idx1] <- colSums( (1/2*repl/v) *
( if(nlregObj$ws$xVar)
v2[,idx1,idx1,drop=FALSE] else 0) +
( (mu*repl-t1)/v)*m2[,idx1,idx1,
drop=FALSE] -
(1/2*si/v^2) *
( if(nlregObj$ws$xVar)
v2[,idx1,idx1,drop=FALSE]
else 0) ) +
crossprod(repl/v*m1[,idx1,drop=FALSE],
m1[,idx1,drop=FALSE]) +
if( nlregObj$ws$xVar)
( crossprod(si/v^3*v1[,idx1,drop=FALSE],
v1[,idx1,drop=FALSE]) -
crossprod(1/2*repl/v^2*v1[,idx1,drop=FALSE],
v1[,idx1,drop=FALSE]) -
crossprod(1/2/v^2*v1[,idx1,drop=FALSE],
sir[,idx1,drop=FALSE]) -
crossprod(1/2/v^2*sir[,idx1,drop=FALSE],
v1[,idx1,drop=FALSE])
)
else 0
if( !is.null(vp) )
{
mat[idx1,idx2] <- if(nlregObj$ws$xVar)
{
( colSums( (1/2*repl/v)*v2[,idx1,idx2,
drop=FALSE] -
(1/2*si/v^2)*v2[,idx1,idx2,
drop=FALSE] ) -
crossprod(1/2*repl/v^2*v1[,idx1,drop=FALSE],
v1[,idx2,drop=FALSE]) +
crossprod(si/v^3*v1[,idx1,drop=FALSE],
v1[,idx2,drop=FALSE]) ) -
crossprod(1/2/v^2*sir[,idx1,drop=FALSE],
v1[,idx2,drop=FALSE] )
}
else if( nlregObj$ws$hom )
- colSums( 1/2/v^2*sir[,idx1,
drop=FALSE]*v1 )
else
- crossprod(1/2/v^2*sir[,idx1,
drop=FALSE],
v1[,,drop=FALSE]) ##
mat[idx2,idx1] <- t(mat[idx1,idx2])
mat[idx2,idx2] <- if(nlregObj$ws$xVar)
( colSums( (repl/2/v)*v2[,idx2,idx2,
drop=FALSE] -
(si/2/v^2)*v2[,idx2,idx2,
drop=FALSE] ) -
crossprod(repl/2/v^2*v1[,idx2,drop=FALSE],
v1[,idx2,drop=FALSE]) +
crossprod(si/v^3*v1[,idx2,drop=FALSE],
v1[,idx2,drop=FALSE])
)
else if(!nlregObj$ws$homVar)
( colSums( (repl/2/v)*v2[,,,drop=FALSE] -
(si/2/v^2)*v2[,,,drop=FALSE] ) -
crossprod(repl/2/v^2*v1[,,drop=FALSE],
v1[,,drop=FALSE]) +
crossprod(si/v^3*v1[,,drop=FALSE],
v1[,,drop=FALSE])
)
else
( sum((repl/2/v) - (si/2/v^2))*v2 -
sum(repl/2/v^2)*v1^2 +
sum(si/v^3)*v1^2
)
}
}
dimnames(mat) <- list(c(names(rc), names(vp)),
c(names(rc), names(vp)))
mat
}
expInfo <- function(object, ...)
UseMethod("expInfo")
expInfo.nlreg <- function(object, par, mu, v, m1=NULL, v1=NULL, ...)
{
nlregObj <- object
rc <- nlregObj$coef
vp <- nlregObj$varPar
of <- nlregObj$offset
repl <- nlregObj$data$repl
if( nlregObj$ws$hoa )
{
md <- nlregObj$ws$md
vd <- nlregObj$ws$vd
}
if( missing(par) ) par <- nlregObj$ws$allPar
if( missing(mu) ) mu <- nlregObj$fitted
if( missing(v) ) v <- nlregObj$weights
.probl <- ( nlregObj$ws$homVar && !is.null(of) )
if(.probl)
.probl <- ( names(of) =="logs" )
# attach(nlregObj$data, warn.conflicts = FALSE) ## 20.05.13
# on.exit( detach(nlregObj$data) )
if( is.null(m1) )
{
tmp <- as.list(par)
tmp <- c(nlregObj$data, tmp)
temp <- if( nlregObj$ws$hoa )
{
formals(md) <- tmp
do.call("md", tmp)
}
else
{
tmp.fun <- Dmean(nlregObj, hessian=FALSE)
formals(tmp.fun) <- tmp
do.call("tmp.fun", tmp)
}
m1 <- attr(temp, "gradient")
m1[!is.finite(m1)] <- 0
}
if(!.probl && is.null(v1))
{
tmp <- as.list(par)
tmp <- c(nlregObj$data, tmp)
temp <- if( nlregObj$ws$hoa )
{
formals(vd) <- tmp
do.call("vd", tmp)
}
else
{
tmp.fun <- Dvar(nlregObj, hessian=FALSE)
formals(tmp.fun) <- tmp
do.call("tmp.fun", tmp)
}
v1 <- attr(temp, "gradient")
v1[!is.finite(v1)] <- 0
}
der <- theta.deriv(nlregObj, par=par, mu=mu, v=v, m1=m1, v1=v1)
if( !nlregObj$ws$missingData )
{
mu <- mu[!duplicated(nlregObj$data$dupl)]
v <- v[!duplicated(nlregObj$data$dupl)]
}
if( nlregObj$ws$homVar ) v <- v[1]
if( .probl )
der2 <- crossprod( der, vecmat(v*repl, der) )
else
{
der2 <- rbind( cbind( diag(v*repl),
if( nlregObj$ws$homVar)
(xx <- as.vector(repl*2*v*mu))
else diag(repl*2*v*mu) ),
cbind( if( nlregObj$ws$homVar) t(xx)
else diag(repl*2*v*mu),
if( nlregObj$ws$homVar )
sum(repl*(2*v^2+4*mu^2*v))
else diag(repl*(2*v^2+4*mu^2*v)) ) )
der2 <- t(der) %*% der2 %*% der
if( nlregObj$ws$homVar || !nlregObj$ws$xVar )
{
der2[1:length(rc),(length(rc)+1):(length(rc)+length(vp))] <- 0
der2[(length(rc)+1):(length(rc)+length(vp)),1:length(rc)] <- 0
}
}
dimnames(der2) <- list(c(names(rc), names(vp)),
c(names(rc), names(vp)))
der2
}
theta.deriv <- function(nlregObj, par, mu, v, m1=NULL, v1=NULL)
{
rc <- nlregObj$coef
vp <- nlregObj$varPar
of <- nlregObj$offset
repl <- nlregObj$data$repl
if( missing(par) ) par <- nlregObj$ws$allPar
if( missing(mu) ) mu <- nlregObj$fitted
if( missing(v) ) v <- nlregObj$weights
.probl <- ( nlregObj$ws$homVar && !is.null(of) )
if(.probl)
.probl <- ( names(of) =="logs" )
# attach(nlregObj$data, warn.conflicts = FALSE) ## 20.05.13
# on.exit( detach(nlregObj$data) )
if( !nlregObj$ws$missingData )
{
mu <- mu[!duplicated(nlregObj$data$dupl)]
v <- v[!duplicated(nlregObj$data$dupl)]
}
if( nlregObj$ws$homVar ) v <- v[1]
if( is.null(m1) )
{
tmp <- as.list(par)
tmp <- c(nlregObj$data, tmp)
temp <- Dmean(nlregObj, hessian=FALSE)
formals(temp) <- tmp
temp <- do.call("temp", tmp)
m1 <- attr(temp, "gradient") ; m1[!is.finite(m1)] <- 0
}
if( !.probl && is.null(v1) )
{
tmp <- as.list(par)
tmp <- c(nlregObj$data, tmp)
temp <- Dvar(nlregObj, hessian=FALSE)
formals(temp) <- tmp
temp <- do.call(temp, tmp)
v1 <- attr(temp, "gradient") ; v1[!is.finite(v1)] <- 0
}
if( .probl )
{
mat <- array(0, dim=c(length(repl), length(rc)))
mat <- 1/v*m1[,,drop=FALSE]
}
else
{
mat <- array(0, dim=c(length(repl)+if( nlregObj$ws$homVar )
1 else length(repl),
length(c(rc,vp))))
idx1 <- 1:length(rc)
if( !is.null(vp) )
idx2 <- length(idx1) + 1:length(vp)
mat[1:length(repl),idx1] <- 1/v*m1[,idx1,drop=FALSE] -
if( nlregObj$ws$xVar )
mu/v^2*v1[,idx1,drop=FALSE] else 0
if( !is.null(vp) )
mat[,idx2] <- rbind(as.matrix(-mu/v^2 * if( nlregObj$ws$xVar )
v1[,idx2,drop=FALSE]
else v1[,,drop=FALSE]),
1/2/v^2* if( nlregObj$ws$xVar )
v1[,idx2,drop=FALSE] else v1[,,drop=FALSE])
if( nlregObj$ws$xVar )
mat[(length(repl)+1):(2*length(repl)), idx1] <-
1/2/v^2*v1[,idx1,drop=FALSE]
}
mat
}
Shat.nlreg <- function(nlregObj1, nlregObj0, par.1, par.0,
mu.1, mu.0, v.1, v.0, m1.1=NULL, m1.0=NULL,
v1.1=NULL, v1.0=NULL)
{
rc <- nlregObj1$coef
vp <- nlregObj1$varPar
of <- nlregObj0$offset
repl <- nlregObj1$data$repl
.probl <- ( nlregObj1$ws$homVar && is.null(nlregObj1$varPar) )
if( missing(par.1) ) par.1 <- nlregObj1$ws$allPar
if( missing(par.0) ) par.0 <- nlregObj0$ws$allPar
if( missing(mu.1) ) mu.1 <- nlregObj1$fitted
if( missing(mu.0) ) mu.0 <- nlregObj0$fitted
if( missing(v.1) ) v.1 <- nlregObj1$weights
if( missing(v.0) ) v.0 <- nlregObj0$weights
# attach(nlregObj1$data, warn.conflicts = FALSE) ## 20.05.13
# on.exit( detach(nlregObj1$data) )
if( is.null(m1.1) || is.null(m1.0) )
{
md <- if( nlregObj1$ws$hoa ) nlregObj1$ws$md
else Dmean(nlregObj1, hessian=FALSE)
tmp <- as.list(nlregObj1$ws$allPar)
tmp <- c(nlregObj1$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.1 <- attr(temp, "gradient") ; m1.1[!is.finite(m1.1)] <- 0
tmp <- as.list(nlregObj0$ws$allPar)
tmp <- c(nlregObj1$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.0 <- attr(temp, "gradient") ; m1.0[!is.finite(m1.0)] <- 0
}
if( is.null(v1.1) || is.null(v1.0) )
{
if( !.probl )
{
vd <- if( nlregObj1$ws$hoa ) nlregObj1$ws$vd
else Dvar(nlregObj1, hessian=FALSE)
tmp <- as.list(nlregObj1$ws$allPar)
tmp <- c(nlregObj1$data, tmp)
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.1 <- attr(temp, "gradient") ; v1.1[!is.finite(v1.1)] <- 0
tmp <- as.list(nlregObj0$ws$allPar)
tmp <- c(nlregObj1$data, tmp)
formals(temp) <- tmp
temp <- do.call("vd", tmp)
v1.0 <- attr(temp, "gradient") ; v1.0[!is.finite(v1.0)] <- 0
}
else v1.0 <- v1.1 <- NULL
}
der.1 <- theta.deriv(nlregObj1, par=par.1, mu=mu.1, v=v.1,
m1=m1.1, v1=v1.1)
der.0 <- theta.deriv(nlregObj1, par=par.0, mu=mu.0, v=v.0,
m1=m1.0, v1=v1.0)
if( !nlregObj1$ws$missingData )
{
mu.1 <- mu.1[!duplicated(nlregObj1$data$dupl)]
mu.0 <- mu.0[!duplicated(nlregObj1$data$dupl)]
v.1 <- v.1[!duplicated(nlregObj1$data$dupl)]
v.0 <- v.0[!duplicated(nlregObj1$data$dupl)]
}
if( nlregObj1$ws$homVar ) { v.1 <- v.1[1] ; v.0 <- v.0[1] }
Sigma <- if( .probl )
diag(v.1*repl)
else
(rbind( cbind( diag(v.1*repl),
if( nlregObj1$ws$homVar)
(xx <- as.vector(repl*2*v.1*mu.1))
else diag(repl*2*v.1*mu.1) ),
cbind( if( nlregObj1$ws$homVar) t(xx)
else diag(repl*2*v.1*mu.1),
if( nlregObj1$ws$homVar )
sum(repl*(2*v.1^2+4*mu.1^2*v.1))
else diag(repl*(2*v.1^2+4*mu.1^2*v.1)) )
))
Shat <- t(der.1) %*% Sigma %*% der.0
if( !.probl )
if( nlregObj1$ws$homVar || !nlregObj1$ws$xVar)
Shat[(length(rc)+1):(length(rc)+length(vp)),1:length(rc)] <-
ifelse( abs(Shat[(length(rc)+1):(length(rc)+length(vp)),
1:length(rc)]) < 1e-10,
0, Shat[(length(rc)+1):(length(rc)+length(vp)),
1:length(rc)] )
dimnames(Shat) <- list(c(names(rc), names(vp)),
c(names(rc), names(vp)))
Shat
}
qhat.nlreg <- function(nlregObj1, nlregObj0, par.1, par.0,
mu.1, mu.0, v.1, v.0, m1.1=NULL, v1.1=NULL)
{
rc <- nlregObj1$coef
vp <- nlregObj1$varPar
of <- nlregObj0$offset
repl <- nlregObj1$data$repl
.probl <- ( nlregObj1$ws$homVar && is.null(nlregObj1$varPar) )
if( missing(par.1) ) par.1 <- nlregObj1$ws$allPar
if( missing(par.0) ) par.0 <- nlregObj0$ws$allPar
if( missing(mu.1) ) mu.1 <- nlregObj1$fitted
if( missing(mu.0) ) mu.0 <- nlregObj0$fitted
if( missing(v.1) ) v.1 <- nlregObj1$weights
if( missing(v.0) ) v.0 <- nlregObj0$weights
# attach(nlregObj1$data, warn.conflicts = FALSE) ## 20.05.13
# on.exit( detach(nlregObj1$data) )
if( is.null(m1.1) )
{
md <- if( nlregObj1$ws$hoa ) nlregObj1$ws$md
else Dmean(nlregObj1, hessian=FALSE)
tmp <- as.list(nlregObj1$ws$allPar)
tmp <- c(nlregObj1$data, tmp)
formals(md) <- tmp
temp <- do.call("md", tmp)
m1.1 <- attr(temp, "gradient") ; m1.1[!is.finite(m1.1)] <- 0
}
if( is.null(v1.1) )
{
if( !.probl )
{
vd <- if( nlregObj1$ws$hoa ) nlregObj1$ws$vd
else Dvar(nlregObj1, hessian=FALSE)
tmp <- as.list(nlregObj1$ws$allPar)
tmp <- c(nlregObj1$data, tmp)
formals(vd) <- tmp
temp <- do.call("vd", tmp)
v1.1 <- attr(temp, "gradient") ; v1.1[!is.finite(v1.1)] <- 0
}
else v1.1 <- NULL
}
der.1 <- theta.deriv(nlregObj1, par=par.1, mu=mu.1, v=v.1,
m1=m1.1, v1=v1.1)
if( !nlregObj1$ws$missingData )
{
mu.1 <- mu.1[!duplicated(nlregObj1$data$dupl)]
mu.0 <- mu.0[!duplicated(nlregObj1$data$dupl)]
v.1 <- v.1[!duplicated(nlregObj1$data$dupl)]
v.0 <- v.0[!duplicated(nlregObj1$data$dupl)]
}
if( nlregObj1$ws$homVar ) { v.1 <- v.1[1] ; v.0 <- v.0[1] }
Sigma <- if( .probl )
diag(v.1*repl)
else
(rbind( cbind( diag(v.1*repl),
if( nlregObj1$ws$homVar)
(xx <- as.vector(repl*2*v.1*mu.1))
else diag(repl*2*v.1*mu.1)),
cbind( if( nlregObj1$ws$homVar) t(xx)
else diag(repl*2*v.1*mu.1),
if( nlregObj1$ws$homVar )
sum(repl*(2*v.1^2+4*mu.1^2*v.1))
else diag(repl*(2*v.1^2+4*mu.1^2*v.1)) )
))
thetah <- mu.1/v.1
thetat <- mu.0/v.0
if( !.probl )
{
thetah <- c(thetah, -1/2/v.1)
thetat <- c(thetat, -1/2/v.0)
}
q <- thetah - thetat
qhat <- t(der.1) %*% Sigma %*% q
if( !.probl )
if( !is.null(nlregObj0$offset) )
{
if( nlregObj1$ws$homVar &&
(names(nlregObj0$offset) == "logs") )
qhat[1:length(rc),1] <-
ifelse( abs(qhat[1:length(rc),1]) < 1e-10,
0, qhat[1:length(rc),1] )
}
dimnames(qhat) <- list(c(names(rc), names(vp)), NULL)
qhat
}
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.