1# These functions are 2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland. 3# All rights reserved. 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30phifun.integrand <- function(charfun, muvec, parlink, earg, 31 cc, xi, ee, ss, sigma2.fuzz = 1e-8) { 32 cnew <- cc / ss 33 34 eta.use <- theta2eta(muvec, link = parlink, earg = earg) # Dangerous 35 chfw <- charfun(x = xi * cnew, eta = eta.use, extra = list()) 36 chfw <- prod(chfw) 37 phi <- chfw * exp(-1i * ee * cnew) * exp(-sigma2.fuzz * (cnew^2) / 2) 38 phi 39} # phifun.integrand 40 41 42 43cdffun.integrand <- 44 function(Object, A.qnorm = 4, wtil, lam, xi, ee, ss, 45 sigma2.fuzz = 1e-8, 46 nodes = sqrt(2) * ghn100, 47 qwts = sqrt(2) * ghw100 / (2 * pi)) { 48 if (!(is(Object, "qrrvglm") || is(Object, "rrvglm"))) 49 stop("argument 'Object' must be a 'qrrvglm' or 'rrvglm' Object") 50 famfun <- Object@family 51 infos <- famfun@infos() 52 charfun <- if (is.logical(infos$charfun) && infos$charfun) 53 famfun@charfun else stop("the family function has no charfun slot") 54 parlink <- linkfun(Object)[1] # All the same, choose the first 55 Earg <- Object@misc$earg[[1]] # Dangerous 56 57 N <- length(nodes) 58 intgrnd <- numeric(N) 59 for (k in 1:N) { 60 curnode <- nodes[k] 61 exptrm <- (exp( 1i * curnode * A.qnorm) - 62 exp(-1i * curnode * wtil)) / (1i * curnode) 63 intgrnd[k] <- phifun.integrand(charfun = charfun, 64 muvec = lam, parlink = parlink, earg = Earg, # Dangerous 65 cc = curnode, xi = xi, 66 ee = ee, ss = ss, sigma2.fuzz = sigma2.fuzz) * 67 exptrm * exp(0.5 * curnode^2) 68 } 69 sum(qwts * intgrnd) 70} # cdffun.integrand 71 72 73 74 75charfun.cqo.cdf <- 76 function( 77 bnu, 78 y0, 79 extra , 80 objfun, 81 Object, 82 Coefs, 83 B1bix, 84 misc.list = list(), 85 Everything, 86 mu.function, 87 88 89 A.qnorm = 4, 90 sigma2.fuzz = 1e-8, lower.latvar = NULL, upper.latvar = NULL, 91 nodes = sqrt(2) * ghn100, 92 qwts = sqrt(2) * ghw100 / (2 * pi)) { 93 wt.mean <- function(x, y, alpha = 0.5) (1 - alpha) * x + alpha * y 94 site.scores <- latvar(Object) 95 range.ss <- range(site.scores) 96 if (is.null(lower.latvar)) 97 lower.latvar <- wt.mean(range.ss[1], range.ss[2], 0.0) 98 if (is.null(upper.latvar)) 99 upper.latvar <- wt.mean(range.ss[1], range.ss[2], 1.0) 100 101 famfun <- Object@family 102 infos <- famfun@infos() 103 charfun <- if (is.logical(infos$charfun) && infos$charfun) 104 famfun@charfun else stop("the family function has no charfun slot") 105 Earg <- Object@misc$earg[[1]] # Dangerous 106 107 uopt <- Opt(Object) 108 tol <- Tol(Object)[1, 1, ] # Interpreted as a variance, not SDs 109 parlink <- linkfun(Object)[1] # All the same, choose the first 110 alf <- theta2eta(Max(Object), parlink, earg = list()) 111 112 113 114 115 nu0 <- bnu 116 qtrm <- (nu0 - uopt) / sqrt(tol) # sqrt(tol), not tol 117 eta <- alf - 0.5 * qtrm^2 118 119 muvec <- eta2theta(eta, parlink, earg = Earg) # Dangerous 120 lam <- muvec 121 122 xi <- -qtrm / sqrt(tol) # sqrt(tol), not tol 123 124 nxi <- sqrt(sum(xi^2)) 125 xi <- xi / nxi 126 127 ee <- sum(lam * xi) 128 varfun <- charfun(eta = eta, # log(muvec), 129 extra = list(), # Object@extra, # For now 130 varfun = TRUE) 131 ss <- sqrt(sum(varfun * xi^2)) # For both poissonff and binomialff 132 w.obs <- sum(y0 * xi) 133 wtil <- (w.obs - ee) / ss 134 if (is.na(wtil)) { 135 prb <- 1 # zz 136 prb <- 0 # zz 137 } else { 138 nrm.prb <- pnorm(wtil) 139 prb <- if (wtil < -A.qnorm) 0 else 140 if ( A.qnorm < wtil) 1 else { 141 prbc <- cdffun.integrand(Object, A.qnorm = A.qnorm, 142 wtil = wtil, lam = lam, xi = xi, 143 ee = ee, ss = ss, 144 sigma2.fuzz = sigma2.fuzz) 145 Re(prbc) 146 } 147 } 148 prb 149} # charfun.cqo.cdf 150 151 152 153 154 155 156 157 158charfun.clo.cdf <- 159 function( 160 bnu, 161 y0, 162 extra , 163 objfun, 164 Object, 165 Coefs, 166 B1bix, 167 misc.list = list(), 168 Everything, 169 mu.function, 170 171 172 A.qnorm = 4, 173 sigma2.fuzz = 1e-8, lower.latvar = NULL, upper.latvar = NULL, 174 nodes = sqrt(2) * ghn100, 175 qwts = sqrt(2) * ghw100 / (2 * pi)) { 176 if (length(bnu) > 1) 177 stop("can only handle rank-1 objects") 178 vfam <- intersect(Object@family@vfamily, c("binomialff", "poissonff")) 179 if (!length(vfam)) 180 stop("only 'poissonff' and 'binomialff' families allowed") 181 all.links <- linkfun(Object) 182 parlink <- all.links[1] # All the same, choose the first 183 canon.link <- switch(vfam, 184 binomialff = all(all.links == "logitlink"), 185 poissonff = all(all.links == "loglink")) 186 if (!canon.link) stop("model does not use the canonical link") # else 187 A.mat <- Coefs@A 188 B1.mat <- Coefs@B1 189 Index.corner <- Object@control$Index.corner # Corner constraints 190 191 192 wt.mean <- function(x, y, alpha = 0.5) (1 - alpha) * x + alpha * y 193 site.scores <- latvar(Object) 194 range.ss <- range(site.scores) 195 if (is.null(lower.latvar)) 196 lower.latvar <- wt.mean(range.ss[1], range.ss[2], 0.0) 197 if (is.null(upper.latvar)) 198 upper.latvar <- wt.mean(range.ss[1], range.ss[2], 1.0) 199 200 famfun <- Object@family 201 infos <- famfun@infos() 202 charfun <- if (is.logical(infos$charfun) && infos$charfun) 203 famfun@charfun else stop("the family function has no charfun slot") 204 Earg <- Object@misc$earg[[1]] # Dangerous 205 206 207 208if (FALSE) { 209 uopt <- Opt(Object) 210 tol <- Tol(Object)[1, 1, ] # Interpreted as a variance, not SDs 211 parlink <- linkfun(Object)[1] # All the same, choose the first 212 alf <- theta2eta(Max(Object), parlink, earg = list()) 213} 214 215 216 217 218 219 220 221 222 223 224 nu0 <- bnu 225 226 227 eta <- Coefs@B1["(Intercept)", ] + A.mat %*% nu0 228 eta <- rbind(c(eta)) # Make sure it is 1 x M 229 230 231 muvec <- eta2theta(eta, parlink, earg = Earg) # Dangerous 232 lam <- muvec 233 234 xi <- A.mat[, 1] 235 236 237 238 nxi <- sqrt(sum(xi^2)) 239 xi <- xi / nxi 240 241 242 243 ee <- sum(lam * xi) 244 varfun <- charfun(eta = eta, # log(muvec), 245 extra = list(), # Object@extra, # For now 246 varfun = TRUE) 247 ss <- sqrt(sum(varfun * xi^2)) # For both poissonff and binomialff 248 w.obs <- sum(y0 * xi) 249 wtil <- (w.obs - ee) / ss 250 if (is.na(wtil)) { 251 prb <- 1 # zz 252 prb <- 0 # zz 253 } else { 254 nrm.prb <- pnorm(wtil) 255 prb <- if (wtil < -A.qnorm) 0 else 256 if ( A.qnorm < wtil) 1 else { 257 prbc <- cdffun.integrand(Object, A.qnorm = A.qnorm, 258 wtil = wtil, lam = lam, xi = xi, 259 ee = ee, ss = ss, 260 sigma2.fuzz = sigma2.fuzz) 261 Re(prbc) 262 } 263 } 264 prb 265} # charfun.clo.cdf 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282fnumat2R <- 283 function(object, 284 refit.model = FALSE) { 285 286 287 288 numat <- latvar(object) # After scaling 289 Rank <- Rank(object) 290 control <- object@control 291 292 M <- npred(object) 293 M1 <- npred(object, type = "one.response") 294 if (M1 > 1) 295 stop("this function only works with M1==1 models") 296 nsmall <- nobs(object) 297 298 X.lm <- model.matrix(object, type = "lm") # Has latvar vars too. 299 etol <- control$eq.tolerances 300 Itol <- control$I.tolerances 301 colx1.index <- control$colx1.index 302 colx2.index <- control$colx2.index 303 NOS <- M / M1 304 305 if (!etol) { 306 index.rc <- iam(NA, NA, M = Rank, both = TRUE) 307 numat2 <- numat[, index.rc$row.index, drop = FALSE] * 308 numat[, index.rc$col.index, drop = FALSE] 309 if (Rank > 1) 310 numat2[, -(1:Rank)] <- numat2[, -(1:Rank)] * 2 311 } # !etol 312 if (etol && !Itol) { 313 numat2 <- numat[, 1:Rank, drop = FALSE] * 314 numat[, 1:Rank, drop = FALSE] # Same as Itol 315 } 316 if (Itol) { 317 numat2 <- numat[, 1:Rank, drop = FALSE] * 318 numat[, 1:Rank, drop = FALSE] 319 } # Itol 320 321 322 323 if (Rank >= 2 && (!etol || (etol && !Itol))) 324 stop("cannot currently handle the given scalings") 325 326 327 328 329 ansA <- kronecker(numat, diag(M)) 330 colnames(ansA) <- param.names("A", NCOL(ansA), skip1 = FALSE) 331 332 333 ansD <- kronecker(numat2, if (Itol || etol) matrix(1, M, 1) else 334 diag(M)) * (-0.5) 335 colnames(ansD) <- param.names("D", NCOL(ansD), skip1 = FALSE) 336 ansx1 <- if (length(colx1.index)) 337 kronecker(X.lm[, colx1.index, drop = FALSE], diag(M)) else 338 NULL # Trivial constraints are assumed. 339 if (length(ansx1)) 340 colnames(ansx1) <- param.names("x1.", NCOL(ansx1), skip1 = FALSE) 341 X.vlm <- cbind(ansA, ansD, ansx1) 342 if (!refit.model) 343 return(X.vlm) 344 345 346 mf <- model.frame(object) 347 mt <- attr(mf, "terms") 348 clist <- vector("list", ncol(X.vlm)) 349 names(clist) <- colnames(X.vlm) 350 for (ii in seq(length(clist))) 351 clist[[ii]] <- diag(M) # Wrong but doesnt matter; trivial constraints 352 353 somejunk <- clist 354 inci <- 0 355 for (ii in seq(length(somejunk))) { 356 inci <- inci + 1 357 somejunk[[ii]] <- inci 358 } 359 attr(numat, "assign") <- somejunk 360 attr(X.vlm, "assign") <- somejunk 361 control$Check.cm.rank <- FALSE 362 control$stepsize <- 1 363 control$half.stepsizing <- FALSE 364 control$Check.rank <- FALSE 365 366 eta.mat.orig <- predict(object) 367 OOO.orig <- object@offset 368 if (!length(OOO.orig) || all(OOO.orig == 0)) 369 OOO.orig <- matrix(0, nsmall, M) 370 371 fit1 <- 372 vglm.fit(x = numat, y = depvar(object), 373 w = c(weights(object, type = "prior")), 374 X.vlm.arg = X.vlm, 375 Xm2 = NULL, Terms = mt, 376 constraints = clist, extra = object@extra, 377 etastart = eta.mat.orig, 378 offset = OOO.orig, family = object@family, 379 control = control) 380 if (fit1$iter > 3) 381 warning("refitted model took an unusually large number of ", 382 "iterations to converge") # Usually 1 iteration is needed 383 return(fit1) 384} # fnumat2R 385 386 387 388 389 390 391forG.calib.qrrvglm <- 392 function(bnu0, # coeff3 = NULL, 393 numerator = c("xi", "eta"), 394 lenb1bix = 1) { 395 396 397 398 numerator <- match.arg(numerator, c("xi", "eta"))[1] 399 400 if (FALSE) { 401 if (!is.null(coeff3)) { 402 if (length(coeff3) != 3 || coeff3[3] > 0) 403 stop("bad input for argument 'coeff3'") 404 beta1 <- coeff3[2] # beta0 <- coeff3[1] is unneeded 405 beta2 <- coeff3[3] 406 } 407 if (is.null(uj)) uj <- -0.5 * beta1 / beta2 408 if (is.null(tolj)) tolj <- 1 / sqrt(-2 * beta2) 409 } # FALSE 410 411 412 413 414 415 Rank <- length(bnu0) # Usually 1, sometimes 2. 416 switch(numerator, 417 "xi" = cbind(rep_len(1, Rank), 2 * bnu0, matrix(0, Rank, lenb1bix)), 418 "eta" = cbind(bnu0, bnu0^2, matrix(1, Rank, lenb1bix))) 419} # forG.calib.qrrvglm 420 421 422 423 424 425 426 427 428dzwald.qrrvglm <- 429 function(bnu0, y0, Object, CoefsObject, 430 B1bix, mu.function) { 431 432 433 434 435 M <- npred(Object) 436 if ((M1 <- npred(Object, type = "one.response")) > 1) 437 stop("this function only works with M1==1 models") 438 nsmall <- nobs(Object) 439 NOS <- M / M1 440 etol <- Object@control$eq.tolerances 441 Itol <- Object@control$I.tolerances 442 Earg <- Object@misc$earg 443 all.links <- linkfun(Object) 444 445 if ((Rank <- Rank(Object)) >= 2 && !Itol) 446 warning("can only handle rank-1 (or rank-2 Itol) objects") 447 448 linkfun <- Object@family@linkfun 449 if (!is.function(linkfun)) 450 stop("could not obtain @linkfun") 451 452 453 vfam <- intersect(Object@family@vfamily, c("binomialff", "poissonff")) 454 if (!length(vfam)) 455 stop("only 'poissonff' and 'binomialff' families allowed") 456 457 canon.link <- switch(vfam, 458 binomialff = all(all.links == "logitlink"), 459 poissonff = all(all.links == "loglink")) 460 if (!canon.link) stop("model does not use the canonical link") # else 461 462 463 464 Omegamat <- vcov(Object) 465 dimn <- colnames(Omegamat) 466 467 468 DD <- 0 # Later becomes a matrix. 469 Gmat <- matrix(0, ncol(Omegamat), Rank) 470 for (spp. in 1:NOS) { 471 index.A.spp. <- spp. + (seq(Rank) - 1) * M # Rank-vector 472 index.D.spp. <- if (Itol || etol) (Rank * M) + seq(Rank) else 473 if (!etol) (Rank * M) + 474 spp. + (seq(Rank * (Rank+1) / 2) - 1) * M # Rank-vector 475 dd <- min(which(substr(dimn, 1, 3) == "x1.")) 476 index.B1.spp. <- dd - 1 + spp. + (seq(1000) - 1) * M 477 index.B1.spp. <- index.B1.spp.[index.B1.spp. <= ncol(Omegamat)] 478 all.index <- c(index.A.spp., index.D.spp., index.B1.spp.) 479 480 481 482 483 484 485 if (!all(should.be.TRUE <- substr(dimn[index.D.spp.], 1, 1) == "D")) 486 stop("encountered a bookkeeping error; failed to pick up the ", 487 "D array coefficients") 488 489 490 491 492 493 uj.jay <- CoefsObject@Optimum[, spp.] # Rank-vector 494 tolj.jay <- if (Rank == 1) sqrt(CoefsObject@Tolerance[, , spp.]) else 495 diag(sqrt(CoefsObject@Tolerance[, , spp.])) # Ditto 496 alphaj.jay <- linkfun(CoefsObject@Maximum, 497 extra = Object@extra)[spp.] # Scalar 498 499 500 eta0.jay <- alphaj.jay - 0.5 * sum(((bnu0 - uj.jay) / tolj.jay)^2) 501 eta0.jay <- rbind(eta0.jay) 502 fv0.jay <- mu.function(eta0.jay, extra = Object@extra) 503 fv0.jay <- c(fv0.jay) # Remove array attributes 504 dTheta.deta0j <- dtheta.deta(fv0.jay, 505 all.links[spp.], 506 earg = Earg[[spp.]]) # Scalar 507 dTheta.deta0j <- c(dTheta.deta0j) # Remove array attributes 508 509 510 xi0.jay <- -(bnu0 - uj.jay) / tolj.jay^2 # Rank-vector 511 DD <- DD + dTheta.deta0j * 512 (cbind(xi0.jay) %*% rbind(xi0.jay)) # More general 513 dxi0.dtheta <- forG.calib.qrrvglm(bnu0, numerator = "xi") # R x dim11 514 deta0.dtheta <- forG.calib.qrrvglm(bnu0, numerator = "eta") # Rxdim11 515 Index.All <- cbind(index.A.spp., index.D.spp., 516 rep_len(index.B1.spp., Rank)) # Just to make sure 517 for (rlocal in seq(Rank)) { 518 Gmat[Index.All[rlocal, ], rlocal] <- 519 Gmat[Index.All[rlocal, ], rlocal] + 520 (y0[1, spp.] - fv0.jay) * dxi0.dtheta[rlocal, ] - 521 xi0.jay[rlocal] * dTheta.deta0j * deta0.dtheta[rlocal, ] 522 } # rlocal 523 524 } # for spp. 525 526 527 528 529 DDinv <- solve(DD) 530 muxf <- diag(Rank) + t(Gmat) %*% Omegamat %*% Gmat %*% DDinv 531 Vmat <- DDinv %*% muxf 532 Vmat 533} # dzwald.qrrvglm 534 535 536 537 538 539 540 541 542 543 544 545 546calibrate.qrrvglm.control <- 547 function(object, 548 trace = FALSE, # passed into optim() 549 method.optim = "BFGS", # passed into optim(method = Method) 550 gridSize = ifelse(Rank == 1, 21, 9), 551 varI.latvar = FALSE, ...) { 552 553 Rank <- object@control$Rank 554 eq.tolerances <- object@control$eq.tolerances 555 if (!is.Numeric(gridSize, positive = TRUE, 556 integer.valued = TRUE, length.arg = 1)) 557 stop("bad input for argument 'gridSize'") 558 if (gridSize < 2) 559 stop("'gridSize' must be >= 2") 560 561 list( 562 trace = as.numeric(trace)[1], 563 method.optim = method.optim, 564 gridSize = gridSize, 565 varI.latvar = as.logical(varI.latvar)[1]) 566} # calibrate.qrrvglm.control 567 568 569 570 571 572 if (!isGeneric("calibrate")) 573 setGeneric("calibrate", 574 function(object, ...) standardGeneric("calibrate")) 575 576 577 578 579 580calibrate.qrrvglm <- 581 function(object, 582 newdata = NULL, 583 type = c("latvar", "predictors", "response", "vcov", 584 "everything"), 585 lr.confint = FALSE, # 20180430 586 cf.confint = FALSE, # 20180602 587 level = 0.95, # 20180430 588 initial.vals = NULL, 589 ...) { 590 591 592 593 594 595 596 se.type <- c("dzwald", "asbefore") # Effectively only the 1st one used 597 Quadratic <- if (is.logical(object@control$Quadratic)) 598 object@control$Quadratic else FALSE # T if CQO, F if CAO 599 600 601 newdata.orig <- newdata 602 if (!length(newdata)) { 603 newdata <- data.frame(depvar(object)) 604 } 605 606 607 if (mode(type) != "character" && mode(type) != "name") 608 type <- as.character(substitute(type)) 609 type <- match.arg(type, c("latvar", "predictors", 610 "response", "vcov", "everything"))[1] 611 612 613 614 615 get.SEs <- Quadratic && type %in% c("vcov", "everything") 616 617 618 619 if (mode(se.type) != "character" && mode(se.type) != "name") 620 se.type <- as.character(substitute(se.type)) 621 se.type <- match.arg(se.type, c("dzwald", "asbefore"))[1] 622 if (se.type == "asbefore") warning("'asbefore' is buggy") 623 624 if (!Quadratic && type == "vcov") 625 stop("cannot have 'type=\"vcov\"' when object is ", 626 "a \"rrvgam\" object") 627 628 629 630 if (!all(weights(object, type = "prior") == 1)) 631 warning("not all the prior weights of 'object' are 1; assuming ", 632 "they all are here") 633 634 635 636 637 638 if (!is.data.frame(newdata)) 639 newdata <- data.frame(newdata) 640 if (!all(object@misc$ynames %in% colnames(newdata))) 641 stop("need the following colnames in 'newdata': ", 642 paste(object@misc$ynames, collapse = ", ")) 643 newdata <- newdata[, object@misc$ynames, drop = FALSE] 644 645 646 647 if (!is.matrix(newdata)) 648 newdata <- as.matrix(newdata) 649 650 651 nn <- nrow(newdata) # Number of sites to calibrate 652 653 654 655 656 obfunct <- slot(object@family, object@misc$criterion) # deviance 657 minimize.obfunct <- 658 if (Quadratic) object@control$min.criterion else 659 TRUE # Logical; TRUE for CAO objects because deviance is minimized 660 if (!is.logical(minimize.obfunct)) 661 stop("object@control$min.criterion is not a logical") 662 optim.control <- calibrate.qrrvglm.control(object = object, ...) 663 664 use.optim.control <- optim.control 665 use.optim.control$method.optim <- 666 use.optim.control$gridSize <- 667 use.optim.control$varI.latvar <- NULL 668 669 670 if ((Rank <- object@control$Rank) > 2) 671 stop("currently can only handle Rank = 1 and Rank = 2") 672 Coefobject <- if (Quadratic) { 673 Coef(object, varI.latvar = optim.control$varI.latvar) 674 } else { 675 Coef(object) 676 } 677 678 679 680 if (!length(initial.vals)) { 681 Lvec <- apply(latvar(object), 2, min) 682 Uvec <- apply(latvar(object), 2, max) 683 initial.vals <- if (Rank == 1) 684 cbind(seq(Lvec, Uvec, length = optim.control$gridSize)) else 685 expand.grid(seq(Lvec[1], Uvec[1], length = optim.control$gridSize), 686 seq(Lvec[2], Uvec[2], length = optim.control$gridSize)) 687 } 688 689 690 691 692 693 694 M <- npred(object) 695 v.simple <- if (Quadratic) { 696 length(object@control$colx1.index) == 1 && 697 names(object@control$colx1.index) == "(Intercept)" && 698 (if (any(names(constraints(object)) == "(Intercept)")) 699 trivial.constraints(constraints(object))["(Intercept)"] == 1 else 700 TRUE) 701 } else { 702 FALSE # To simplify things for "rrvgam" objects 703 } 704 B1bix <- if (v.simple) { 705 matrix(Coefobject@B1, nn, M, byrow = TRUE) 706 } else { 707 Xlm <- predict.vlm(as(object, "vglm"), # object, 708 newdata = newdata.orig, 709 type = "Xlm") 710 if (se.type == "dzwald" && (type == "everything" || type == "vcov")) 711 stop("only noRRR = ~ 1 models are handled for ", 712 "type = 'everything' or type = 'vcov'") 713 714 Xlm[, names(object@control$colx1.index), drop = FALSE] %*% 715 (if (Quadratic) Coefobject@B1 else object@coefficients[1:M]) 716 } # !v.simple 717 718 719 720 721 722 723 724 objfun1 <- function(lv1val, 725 x = NULL, y, w = 1, extraargs) { 726 ans <- sum(c(w) * extraargs$Obfunction( 727 bnu = c(lv1val), 728 y0 = y, 729 extra = extraargs$object.extra, 730 objfun = extraargs$obfunct, 731 Object = extraargs$Object, # Needed for "rrvgam" objects 732 Coefs = extraargs$Coefobject, 733 B1bix = extraargs$B1bix, 734 misc.list = extraargs$object.misc, 735 Everything = FALSE, 736 mu.function = extraargs$mu.function)) 737 ans 738 } # objfun1 739 740 objfun2 <- function(lv1val, lv2val, 741 x = NULL, y, w = 1, extraargs) { 742 ans <- sum(c(w) * extraargs$Obfunction( 743 bnu = c(lv1val, lv2val), 744 y0 = y, 745 extra = extraargs$object.extra, 746 objfun = extraargs$obfunct, 747 Object = extraargs$Object, # Needed for "rrvgam" objects 748 Coefs = extraargs$Coefobject, 749 B1bix = extraargs$B1bix, 750 misc.list = extraargs$object.misc, 751 Everything = FALSE, 752 mu.function = extraargs$mu.function)) 753 ans 754 } # objfun2 755 756 757 choose.fun <- if (Quadratic) .my.calib.objfunction.qrrvglm else 758 .my.calib.objfunction.rrvgam 759 mu.function <- slot(object@family, "linkinv") 760 wvec <- 1 # zz; Assumed here 761 mylist <- 762 list(object.extra = object@extra, 763 Obfunction = choose.fun, # e.g. .my.calib.objfunction.qrrvglm 764 Coefobject = Coefobject, 765 B1bix = NA, # Will be replaced below 766 obfunct = obfunct, # deviance 767 object.misc = object@misc, 768 Object = if (Quadratic) 777 else object, 769 mu.function = mu.function) 770 771 772 773 init.vals <- matrix(NA_real_, nn, Rank) 774 for (i1 in 1:nn) { 775 if (optim.control$trace) 776 cat("Grid searching initial values for observation", 777 i1, "-----------------\n") 778 779 780 y0 <- newdata[i1, , drop = FALSE] # drop added 20150624 781 mylist$B1bix <- B1bix[i1, ] 782 try.this <- if (Rank == 1) 783 grid.search(initial.vals[, 1], 784 objfun = objfun1, y = y0 , w = wvec, 785 ret.objfun = TRUE, 786 maximize = !minimize.obfunct, # Most general. 787 extraargs = mylist) else 788 grid.search2(initial.vals[, 1], initial.vals[, 2], 789 objfun = objfun2, y = y0, w = wvec, 790 ret.objfun = TRUE, 791 maximize = !minimize.obfunct, # Most general. 792 extraargs = mylist) 793 lv1.init <- try.this[if (Rank == 1) "Value" else "Value1"] 794 lv2.init <- if (Rank >= 2) try.this["Value2"] else NULL 795 796 init.vals[i1, ] <- c(lv1.init, lv2.init) 797 } # for i1 798 799 800 801 802 803 804 BestOFpar <- matrix(NA_real_, nn, Rank) 805 BestOFvalues <- rep(NA_real_, nn) # Best OF objective function values 806 807 for (i1 in 1:nn) { 808 if (optim.control$trace) { 809 cat("\nOptimizing for observation", i1, "-----------------\n") 810 flush.console() 811 } 812 813 ans <- 814 optim(par = init.vals[i1, ], 815 fn = choose.fun, 816 method = optim.control$method.optim, # "BFGS" or "CG" or... 817 control = c(fnscale = ifelse(minimize.obfunct, 1, -1), 818 use.optim.control), # as.vector() needed 819 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 820 extra = object@extra, 821 objfun = obfunct, # deviance 822 Object = if (Quadratic) 777 else object, 823 Coefs = Coefobject, 824 B1bix = B1bix[i1, , drop = FALSE], 825 misc.list = object@misc, 826 Everything = FALSE, # Differs from Step 6. below 827 mu.function = mu.function) 828 829 if (optim.control$trace) { 830 if (ans$convergence == 0) 831 cat("Successful convergence\n") else 832 cat("Unsuccessful convergence\n") 833 flush.console() 834 } 835 if (ans$convergence == 0) { 836 BestOFpar[i1, ] <- ans$par 837 BestOFvalues[i1] <- ans$value 838 } # else do nothing since NA_real_ signifies convergence failure 839 } # for i1 840 841 842 843 844 845 846 847 prettyCQO <- function(BestOFpar, newdata, Rank) { 848 if (Rank == 1) { 849 BestOFpar <- c(BestOFpar) # Drop the dimension 850 if (!is.null(dimnames(newdata)[[1]])) { 851 names(BestOFpar) <- dimnames(newdata)[[1]] 852 } 853 } else { 854 dimnames(BestOFpar) <- list(dimnames(newdata)[[1]], 855 param.names("latvar", Rank, skip1 = TRUE)) 856 } 857 BestOFpar 858 } # prettyCQO 859 860 BestOFpar <- prettyCQO(BestOFpar, newdata, Rank) 861 attr(BestOFpar,"objectiveFunction") <- 862 prettyCQO(BestOFvalues, newdata, Rank = 1) 863 if (type == "latvar" && (!cf.confint && !lr.confint)) { 864 return(BestOFpar) 865 } 866 867 868 869 870 871 872 873 if (lr.confint && Rank > 1) { 874 warning("argument 'lr.confint' should only be TRUE if Rank == 1. ", 875 "Setting 'lr.confint = FALSE'.") 876 lr.confint <- FALSE 877 } 878 879 880 if (lr.confint && !(type %in% c("latvar", "everything"))) { 881 warning("argument 'lr.confint' should only be TRUE if ", 882 "'type = \"latvar\"' or 'type = \"everything\"'. ", 883 "Setting 'lr.confint = FALSE'.") 884 lr.confint <- FALSE 885 } 886 if (lr.confint && Rank == 1) { 887 888 format.perc <- function(probs, digits) 889 paste(format(100 * probs, trim = TRUE, scientific = FALSE, 890 digits = digits), "%") 891 aa <- (1 - level) / 2 892 aa <- c(aa, 1 - aa) 893 894 pct <- format.perc(aa, 3) 895 cimat <- array(NA, dim = c(nn, 2L), 896 dimnames = list(dimnames(newdata)[[1]], pct)) 897 898 for (i1 in 1:nn) { 899 if (optim.control$trace) { 900 cat("\nSolving for the roots for obsn", i1, "---------------\n") 901 flush.console() 902 } 903 904 905 foo1.lhs.rhs <- 906 function(bnu, 907 y0, extra = NULL, 908 objfun, Coefs, 909 Object, # Not needed 910 B1bix, 911 misc.list, 912 Everything = FALSE, 913 mu.function, 914 BestOFvalues = NA, 915 level = 0.95, 916 criterion.arg = "loglikelihood") { 917 if (!(criterion.arg %in% c("deviance", "loglikelihood"))) 918 stop("'criterion.arg' must be 'deviance' or 'loglikelihood'") 919 ifelse(criterion.arg == "deviance", 1, -2) * 920 (-BestOFvalues + 921 .my.calib.objfunction.qrrvglm(bnu = bnu, 922 y0 = y0, 923 extra = extra, 924 objfun = objfun, 925 Object = Object, 926 Coefs = Coefs, 927 B1bix = B1bix, 928 misc.list = misc.list, 929 Everything = Everything, 930 mu.function = mu.function)) - 931 qchisq(level, df = 1) 932 } 933 934 935 for (Side in 1:2) { 936 ans.lhs.rhs <- 937 uniroot(f = foo1.lhs.rhs, 938 interval = if (Side == 1) c(Lvec, BestOFpar[i1]) else 939 c(BestOFpar[i1], Uvec), 940 extendInt = ifelse(Side == 1, "downX", "upX"), 941 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 942 extra = object@extra, 943 objfun = obfunct, 944 Coefs = Coefobject, 945 Object = 777, # object, 946 B1bix = B1bix[i1, , drop = FALSE], 947 misc.list = object@misc, 948 Everything = FALSE, 949 mu.function = mu.function, 950 BestOFvalues = BestOFvalues[i1], level = level, 951 criterion.arg = object@misc$criterion) 952 cimat[i1, Side] <- ans.lhs.rhs$root 953 } # Side 954 } # for i1 955 956 if (type == "latvar") 957 return(cbind(estimate = BestOFpar, 958 objfun = BestOFvalues, 959 cimat)) 960 } # if (lr.confint && Rank == 1) 961 962 963 964 965 966 967 968 969 970 if (cf.confint && Rank > 1) { 971 warning("argument 'cf.confint' should only be TRUE if Rank == 1. ", 972 "Setting 'cf.confint = FALSE'.") 973 cf.confint <- FALSE 974 } 975 976 977 if (cf.confint && !(type %in% c("latvar", "everything"))) { 978 warning("argument 'cf.confint' should only be TRUE if ", 979 "'type = \"latvar\"' or 'type = \"everything\"'. ", 980 "Setting 'cf.confint = FALSE'.") 981 cf.confint <- FALSE 982 } 983 if (cf.confint && Rank == 1) { 984 985 format.perc <- function(probs, digits) 986 paste(format(100 * probs, trim = TRUE, scientific = FALSE, 987 digits = digits), "%") 988 aa <- (1 - level) / 2 989 aa <- c(aa, 1 - aa) 990 pct <- format.perc(aa, 3) 991 cimat2 <- array(NA, dim = c(nn, 2L), 992 dimnames = list(dimnames(newdata)[[1]], pct)) 993 994 995 996 for (i1 in 1:nn) { 997 if (optim.control$trace) { 998 cat("\nSolving for the roots for obsn", i1, "---------------\n") 999 flush.console() 1000 } 1001 1002 1003 foo2.lhs.rhs <- 1004 function(bnu, 1005 y0, extra = NULL, 1006 objfun, Coefs, 1007 Object, # Not needed 1008 B1bix, 1009 misc.list, 1010 Everything = FALSE, 1011 mu.function, 1012 BestOFvalues = NA, 1013 pr.level = c(0.05, 0.95)[1] 1014 ) { 1015 charfun.cqo.cdf(bnu = bnu, 1016 y0 = y0, 1017 extra = extra, 1018 objfun = objfun, 1019 Coefs = Coefs, 1020 Object = Object, 1021 B1bix = B1bix, 1022 misc.list = misc.list, 1023 Everything = Everything, 1024 mu.function = mu.function 1025 ) - 1026 pr.level 1027 } 1028 1029 1030 1031 1032 1033 1034 1035 for (Side in 1:2) { 1036 ans.lhs.rhs <- 1037 uniroot(f = foo2.lhs.rhs, 1038 interval = if (Side == 1) c(Lvec, BestOFpar[i1]) else 1039 c(BestOFpar[i1], Uvec), 1040 extendInt = "yes", # Might be worse than above. 1041 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 1042 extra = object@extra, 1043 objfun = obfunct, 1044 Coefs = Coefobject, 1045 Object = object, 1046 B1bix = B1bix[i1, , drop = FALSE], 1047 misc.list = object@misc, 1048 Everything = FALSE, 1049 mu.function = mu.function, 1050 BestOFvalues = BestOFvalues[i1], 1051 pr.level = ifelse(Side == 1, aa[1], aa[2]) 1052 ) 1053 cimat2[i1, Side] <- ans.lhs.rhs$root 1054 } # Side 1055 } # for i1 1056 1057 vecTF <- cimat2[, 2] < cimat2[, 1] 1058 if (any(vecTF)) { 1059 temp <- cimat2[vecTF, 1] 1060 cimat2[vecTF, 1] <- cimat2[vecTF, 2] 1061 cimat2[vecTF, 2] <- temp 1062 } 1063 if (type == "latvar") 1064 return(cbind(estimate = BestOFpar, 1065 objfun = BestOFvalues, 1066 cimat2)) 1067 1068 } # if (cf.confint && Rank == 1) 1069 1070 1071 1072 1073 1074 1075 1076 1077 etaValues <- matrix(NA_real_, nn, M) 1078 muValues <- matrix(NA_real_, nn, ncol(fitted(object))) 1079 vcValues <- if (get.SEs) array(0, c(Rank, Rank, nn)) else NULL 1080 1081 1082 if (optim.control$trace) 1083 cat("\n") 1084 1085 for (i1 in 1:nn) { 1086 if (optim.control$trace) { 1087 cat("Evaluating quantities for observation", i1, 1088 "-----------------\n") 1089 flush.console() 1090 } 1091 ans5 <- choose.fun( 1092 bnu = if (Rank == 1) BestOFpar[i1] else BestOFpar[i1, ], 1093 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 1094 extra = object@extra, 1095 objfun = obfunct, # deviance 1096 Object = if (Quadratic) 777 else object, 1097 Coefs = Coefobject, 1098 B1bix = B1bix[i1, , drop = FALSE], 1099 misc.list = object@misc, 1100 Everything = TRUE, # Differs from Step 3. 1101 mu.function = mu.function) 1102 1103 muValues[i1, ] <- ans5$mu 1104 etaValues[i1, ] <- ans5$eta 1105 1106 1107 if (get.SEs) { 1108 vcValues[, , i1] <- if (se.type == "dzwald") 1109 dzwald.qrrvglm( 1110 bnu0 = if (Rank == 1) BestOFpar[i1] else BestOFpar[i1, ], 1111 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 1112 Object = object, 1113 CoefsObject = Coefobject, 1114 B1bix = B1bix[i1, , drop = FALSE], 1115 mu.function = mu.function) else 1116 ans5$vcmat # Might be NULL, e.g., "rrvgam" 1117 } # if (get.SEs) 1118 } # for i1 1119 1120 1121 1122 1123 1124 1125 1126 dimnames(muValues) <- dimnames(newdata) 1127 dimnames(etaValues) <- list(dimnames(newdata)[[1]], 1128 dimnames(object@predictors)[[2]]) 1129 if (get.SEs) 1130 dimnames(vcValues) <- list(as.character(1:Rank), 1131 as.character(1:Rank), 1132 dimnames(newdata)[[1]]) 1133 1134 switch(type, 1135 latvar = BestOFpar, # Done already, so not really needed 1136 predictors = etaValues, 1137 response = muValues, 1138 vcov = vcValues, 1139 everything = list( 1140 latvar = BestOFpar, 1141 predictors = etaValues, 1142 response = muValues, 1143 vcov = vcValues, 1144 lr.confint = if (lr.confint) 1145 cbind(estimate = BestOFpar, 1146 objfun = BestOFvalues, 1147 cimat) else NULL, 1148 cf.confint = if (cf.confint) 1149 cbind(estimate = BestOFpar, 1150 objfun = BestOFvalues, 1151 cimat2) else NULL) 1152 ) 1153} # calibrate.qrrvglm 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165.my.calib.objfunction.qrrvglm <- 1166 function(bnu, y0, extra = NULL, 1167 objfun, Coefs, Object, 1168 B1bix, 1169 misc.list, 1170 Everything = TRUE, 1171 mu.function) { 1172 1173 1174 bnumat <- cbind(bnu) 1175 Rank <- length(bnu) 1176 eta <- cbind(c(B1bix)) + Coefs@A %*% bnumat 1177 M <- misc.list$M 1178 check.eta <- matrix(0, M, 1) 1179 for (ss in 1:M) { 1180 temp <- Coefs@D[, , ss, drop = FALSE] 1181 dim(temp) <- dim(temp)[1:2] # c(Rank, Rank) 1182 eta[ss, 1] <- eta[ss, 1] + t(bnumat) %*% temp %*% bnumat 1183 1184 if (FALSE) { 1185 warning("this line is wrong:") 1186 alf <- loglink(Coefs@Maximum[ss]) # zz get the link function 1187 tolmat <- Coefs@Tolerance[, , ss, drop = FALSE] 1188 check.eta[ss, 1] <- alf - 0.5 * t(bnumat) %*% 1189 solve(tolmat) %*% bnumat 1190 } # FALSE 1191 } # for ss 1192 eta <- matrix(eta, 1, M, byrow = TRUE) 1193 mu <- matrix(mu.function(eta, extra = extra), nrow = 1) 1194 obvalue <- objfun(mu = mu, y = y0, 1195 w = 1, # ignore prior.weights on the object 1196 residuals = FALSE, eta = eta, extra = extra) 1197 if (Everything) { 1198 vcmat <- diag(Rank) 1199 if (FALSE && M == NCOL(mu)) { 1200 for (ss in 1:M) { 1201 vec1 <- cbind(Coefs@A[ss, ]) + 1202 2 * matrix(Coefs@D[, , ss], Rank, Rank) %*% bnumat 1203 vcmat <- vcmat + mu[1, ss] * vec1 %*% t(vec1) 1204 } # ss 1205 } # if (M == NCOL(mu)) 1206 vcmat <- solve(vcmat) 1207 } else { 1208 vcmat <- NULL 1209 } 1210 if (Everything) 1211 list(eta = eta, 1212 mu = mu, 1213 obvalue = obvalue, 1214 vcmat = vcmat) else 1215 obvalue 1216} # .my.calib.objfunction.qrrvglm 1217 1218 1219 1220 1221 1222 1223.my.calib.objfunction.rrvgam <- 1224 function(bnu, y0, extra = NULL, 1225 objfun, 1226 Object, # Needed for "rrvgam" objects 1227 Coefs, 1228 B1bix, # Actually not needed here 1229 misc.list, 1230 Everything = TRUE, 1231 mu.function) { 1232 Rank <- length(bnu) 1233 NOS <- Coefs@NOS 1234 eta <- matrix(NA_real_, 1, NOS) 1235 for (jlocal in 1:NOS) { 1236 eta[1, jlocal] <- predictrrvgam(Object, grid = bnu, 1237 sppno = jlocal, Rank = Rank, deriv = 0)$yvals 1238 } 1239 mu <- rbind(mu.function(eta, extra)) # Make sure it has one row 1240 obvalue <- objfun(mu = mu, y = y0, 1241 w = 1, # ignore prior.weights on the object 1242 residuals = FALSE, eta = eta, extra = extra) 1243 vcmat <- NULL # No theory as of yet to compute the vcmat 1244 if (Everything) 1245 list(eta = eta, 1246 mu = mu, 1247 obvalue = obvalue, 1248 vcmat = vcmat) else 1249 obvalue 1250} # .my.calib.objfunction.rrvgam 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264forG.calib.rrvglm <- 1265 function(bnu0, 1266 numerator = c("xi", "eta"), 1267 lenb1bix = 1) { 1268 1269 numerator <- match.arg(numerator, c("xi", "eta"))[1] 1270 1271 Rank <- length(bnu0) # Usually 1, sometimes 2. 1272 switch(numerator, 1273 "xi" = cbind(rep_len(1, Rank), matrix(0, Rank, lenb1bix)), 1274 "eta" = cbind(bnu0, matrix(1, Rank, lenb1bix))) 1275} # forG.calib.rrvglm 1276 1277 1278 1279 1280 1281 1282dzwald.rrvglm <- 1283 function(bnu0, y0, # extra = NULL, objfun, 1284 Object, CoefsObject, 1285 B1bix, 1286 mu.function 1287 ) { 1288 1289 1290 M <- npred(Object) 1291 if ((M1 <- npred(Object, type = "one.response")) > 1) 1292 stop("this function only works with M1==1 models") 1293 NOS <- M / M1 1294 Earg <- Object@misc$earg 1295 all.links <- linkfun(Object) 1296 1297 1298 1299 1300 Rank <- Rank(Object) 1301 1302 linkfun <- Object@family@linkfun 1303 if (!is.function(linkfun)) 1304 stop("could not obtain @linkfun") 1305 1306 1307 vfam <- intersect(Object@family@vfamily, c("binomialff", "poissonff")) 1308 if (!length(vfam)) 1309 stop("only 'poissonff' and 'binomialff' families allowed") 1310 1311 canon.link <- switch(vfam, 1312 binomialff = all(all.links == "logitlink"), 1313 poissonff = all(all.links == "loglink")) 1314 if (!canon.link) stop("model does not use the canonical link") # else 1315 1316 1317 1318 Omegamat <- vcov(Object) # Numerical problems might occur to get this. 1319 dimn <- colnames(Omegamat) 1320 A.mat <- CoefsObject@A 1321 B1.mat <- CoefsObject@B1 1322 Index.corner <- Object@control$Index.corner # Corner constraints 1323 1324 DD <- 0 # Later becomes a matrix. 1325 Gmat <- matrix(0, ncol(Omegamat), Rank) 1326 icounter <- 0 # Number of rows of \bigtilde{\bA}. 1327 for (spp. in 1:NOS) { 1328 index.A.spp. <- if (any(spp. == Index.corner)) NULL else { 1329 icounter <- icounter + 1 1330 icounter + (M - Rank) * (seq(Rank) - 1) 1331 } 1332 index.D.spp. <- NULL 1333 dd <- max(which(substr(dimn, 1, 8) == "I(latvar")) 1334 index.B1.spp. <- dd + spp. + (seq(nrow(B1.mat)) - 1) * M 1335 1336 1337 all.index <- c(index.A.spp., index.D.spp., index.B1.spp.) 1338 1339 1340 1341 1342 1343 1344 1345 alphaj.jay <- CoefsObject@B1["(Intercept)", spp.] # Scalar 1346 1347 1348 eta0.jay <- alphaj.jay + A.mat[spp., , drop = FALSE] %*% bnu0 1349 eta0.jay <- rbind(eta0.jay) 1350 fv0.jay <- mu.function(eta0.jay, extra = Object@extra) 1351 fv0.jay <- c(fv0.jay) # Remove array attributes 1352 dTheta.deta0j <- dtheta.deta(fv0.jay, 1353 all.links[spp.], 1354 earg = Earg[[spp.]]) # Scalar 1355 dTheta.deta0j <- c(dTheta.deta0j) # Remove array attributes 1356 1357 1358 xi0.jay <- A.mat[spp., ] # Rank-vector 1359 DD <- DD + dTheta.deta0j * 1360 (cbind(xi0.jay) %*% rbind(xi0.jay)) # More general 1361 dxi0.dtheta <- forG.calib.rrvglm(bnu0, numerator = "xi") # R x dim11 1362 deta0.dtheta <- forG.calib.rrvglm(bnu0, numerator = "eta") # Rxdim11 1363 Index.All <- cbind(index.A.spp., index.D.spp., 1364 rep_len(index.B1.spp., Rank)) # Just to make sure 1365 1366 if (!is.null(index.A.spp.)) 1367 for (rlocal in seq(Rank)) { 1368 Gmat[Index.All[rlocal, ], rlocal] <- 1369 Gmat[Index.All[rlocal, ], rlocal] + 1370 (y0[1, spp.] - fv0.jay) * dxi0.dtheta[rlocal, ] - 1371 xi0.jay[rlocal] * dTheta.deta0j * deta0.dtheta[rlocal, ] 1372 } # rlocal 1373 1374 } # for spp. 1375 1376 1377 1378 1379 DDinv <- solve(DD) 1380 muxf <- diag(Rank) + t(Gmat) %*% Omegamat %*% Gmat %*% DDinv 1381 Vmat <- DDinv %*% muxf 1382 Vmat 1383} # dzwald.rrvglm 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393calibrate.rrvglm <- 1394 function(object, 1395 newdata = NULL, 1396 type = c("latvar", "predictors", "response", "vcov", 1397 "everything"), 1398 lr.confint = FALSE, # 20180427 1399 cf.confint = FALSE, # 20180604 1400 level = 0.95, # 20180428 1401 initial.vals = NULL, # For one observation only 1402 ...) { 1403 se.type <- c("dzwald", "asbefore") # Effectively only the 1st one used 1404 1405 1406 Quadratic <- FALSE # Because this function was adapted from CQO code. 1407 newdata.orig <- newdata 1408 if (!length(newdata)) { 1409 newdata <- data.frame(depvar(object)) 1410 } 1411 1412 if (mode(type) != "character" && mode(type) != "name") 1413 type <- as.character(substitute(type)) 1414 type <- match.arg(type, c("latvar", "predictors", 1415 "response", "vcov", "everything"))[1] 1416 get.SEs <- type %in% c("vcov", "everything") 1417 if (mode(se.type) != "character" && mode(se.type) != "name") 1418 se.type <- as.character(substitute(se.type)) 1419 se.type <- match.arg(se.type, c("dzwald", "asbefore"))[1] 1420 1421 if (!all(weights(object, type = "prior") == 1)) 1422 warning("not all the prior weights of 'object' are 1; assuming ", 1423 "they all are here") 1424 1425 1426 1427 1428 if (!is.data.frame(newdata)) 1429 newdata <- data.frame(newdata) 1430 if (!all(object@misc$ynames %in% colnames(newdata))) 1431 stop("need the following colnames in 'newdata': ", 1432 paste(object@misc$ynames, collapse = ", ")) 1433 newdata <- newdata[, object@misc$ynames, drop = FALSE] 1434 1435 1436 if (!is.matrix(newdata)) 1437 newdata <- as.matrix(newdata) 1438 1439 1440 nn <- nrow(newdata) # Number of sites to calibrate 1441 1442 1443 obfunct <- slot(object@family, object@misc$criterion) 1444 minimize.obfunct <- object@control$min.criterion # deviance 1445 if (!is.logical(minimize.obfunct)) 1446 stop("object@control$min.criterion is not a logical") 1447 minimize.obfunct <- as.vector(minimize.obfunct) 1448 optim.control <- calibrate.rrvglm.control(object = object, ...) 1449 1450 use.optim.control <- optim.control 1451 use.optim.control$method.optim <- 1452 use.optim.control$gridSize <- 1453 use.optim.control$varI.latvar <- NULL 1454 1455 1456 if ((Rank <- object@control$Rank) > 3) 1457 stop("currently can only handle Rank = 1, 2 and 3") 1458 Coefobject <- if (Quadratic) { 1459 Coef(object, varI.latvar = optim.control$varI.latvar) 1460 } else { 1461 Coef(object) 1462 } 1463 1464 1465 1466 1467 1468 if (!length(initial.vals)) { 1469 Lvec <- apply(latvar(object), 2, min) 1470 Uvec <- apply(latvar(object), 2, max) 1471 initial.vals <- if (Rank == 1) 1472 cbind(seq(Lvec, Uvec, length = optim.control$gridSize)) else 1473 if (Rank == 2) 1474 expand.grid(seq(Lvec[1], Uvec[1], length = optim.control$gridSize), 1475 seq(Lvec[2], Uvec[2], length = optim.control$gridSize)) else 1476 expand.grid(seq(Lvec[1], Uvec[1], length = optim.control$gridSize), 1477 seq(Lvec[2], Uvec[2], length = optim.control$gridSize), 1478 seq(Lvec[3], Uvec[3], length = optim.control$gridSize)) 1479 } # !length(initial.vals) 1480 1481 1482 1483 1484 M <- npred(object) 1485 v.simple <- length(object@control$colx1.index) == 1 && 1486 names(object@control$colx1.index) == "(Intercept)" && 1487 (if (any(names(constraints(object)) == "(Intercept)")) 1488 trivial.constraints(constraints(object))["(Intercept)"] == 1 else 1489 TRUE) 1490 B1bix <- if (v.simple) { 1491 matrix(Coefobject@B1, nn, M, byrow = TRUE) 1492 } else { 1493 Xlm <- predict.vlm(as(object, "vglm"), # object, 1494 newdata = newdata.orig, 1495 type = "Xlm") 1496 if (NROW(Xlm) != nn) 1497 warning("NROW(Xlm) and ", nn, " are unequal") 1498 1499 if (se.type == "dzwald" && (type == "everything" || type == "vcov")) 1500 stop("only noRRR = ~ 1 models are handled for ", 1501 "type = 'everything' or type = 'vcov'") 1502 1503 Xlm[, names(object@control$colx1.index), drop = FALSE] %*% 1504 Coefobject@B1 1505 } # !v.simple 1506 1507 1508 1509 1510 1511 objfun1 <- function(lv1val, 1512 x = NULL, y, w = 1, extraargs) { 1513 ans <- sum(c(w) * extraargs$Obfunction( 1514 bnu = c(lv1val), 1515 y0 = y, 1516 extra = extraargs$object.extra, 1517 objfun = extraargs$obfunct, 1518 Object = extraargs$Object, 1519 Coefs = extraargs$Coefobject, 1520 B1bix = extraargs$B1bix, 1521 misc.list = extraargs$object.misc, 1522 Everything = FALSE, 1523 mu.function = extraargs$mu.function)) 1524 ans 1525 } 1526 1527 objfun2 <- function(lv1val, lv2val, 1528 x = NULL, y, w = 1, extraargs) { 1529 ans <- sum(c(w) * extraargs$Obfunction( 1530 bnu = c(lv1val, lv2val), 1531 y0 = y, 1532 extra = extraargs$object.extra, 1533 objfun = extraargs$obfunct, 1534 Object = extraargs$Object, 1535 Coefs = extraargs$Coefobject, 1536 B1bix = extraargs$B1bix, 1537 misc.list = extraargs$object.misc, 1538 Everything = FALSE, 1539 mu.function = extraargs$mu.function)) 1540 ans 1541 } 1542 1543 objfun3 <- function(lv1val, lv2val, lv3val, 1544 x = NULL, y, w = 1, extraargs) { 1545 ans <- sum(c(w) * extraargs$Obfunction( 1546 bnu = c(lv1val, lv2val, lv3val), 1547 y0 = y, 1548 extra = extraargs$object.extra, 1549 objfun = extraargs$obfunct, 1550 Object = extraargs$Object, 1551 Coefs = extraargs$Coefobject, 1552 B1bix = extraargs$B1bix, 1553 misc.list = extraargs$object.misc, 1554 Everything = FALSE, 1555 mu.function = extraargs$mu.function)) 1556 ans 1557 } 1558 1559 1560 mu.function <- slot(object@family, "linkinv") 1561 wvec <- 1 # zz; Assumed here 1562 mylist <- 1563 list(object.extra = object@extra, 1564 Obfunction = .my.calib.objfunction.rrvglm, 1565 Coefobject = Coefobject, 1566 B1bix = NA, # Will be replaced below 1567 obfunct = obfunct, 1568 object.misc = object@misc, 1569 Object = 777, # object, 1570 mu.function = mu.function) 1571 1572 1573 init.vals <- matrix(NA_real_, nn, Rank) 1574 for (i1 in 1:nn) { 1575 if (optim.control$trace) 1576 cat("Grid searching initial values for observation", 1577 i1, "-----------------\n") 1578 1579 1580 y0 <- newdata[i1, , drop = FALSE] # drop added 20150624 1581 mylist$B1bix <- B1bix[i1, ] 1582 try.this <- if (Rank == 1) 1583 grid.search(initial.vals[, 1], 1584 objfun = objfun1, y = y0 , w = wvec, 1585 ret.objfun = TRUE, 1586 maximize = !minimize.obfunct, # Most general. 1587 extraargs = mylist) else if (Rank == 2) 1588 grid.search2(initial.vals[, 1], initial.vals[, 2], 1589 objfun = objfun2, y = y0, w = wvec, 1590 ret.objfun = TRUE, 1591 maximize = !minimize.obfunct, # Most general. 1592 extraargs = mylist) else 1593 grid.search3(initial.vals[, 1], initial.vals[, 2], 1594 initial.vals[, 3], 1595 objfun = objfun3, y = y0, w = wvec, 1596 ret.objfun = TRUE, 1597 maximize = !minimize.obfunct, # Most general. 1598 extraargs = mylist) 1599 lv1.init <- try.this[if (Rank == 1) "Value" else "Value1"] 1600 lv2.init <- if (Rank >= 2) try.this["Value2"] else NULL 1601 lv3.init <- if (Rank >= 3) try.this["Value3"] else NULL 1602 1603 init.vals[i1, ] <- c(lv1.init, lv2.init, lv3.init) 1604 } # for i1 1605 1606 1607 1608 1609 1610 1611 BestOFpar <- matrix(NA_real_, nn, Rank) 1612 BestOFvalues <- rep(NA_real_, nn) # Best OF objective function values 1613 1614 for (i1 in 1:nn) { 1615 if (optim.control$trace) { 1616 cat("\nOptimizing for observation", i1, "-----------------\n") 1617 flush.console() 1618 } 1619 1620 ans <- 1621 optim(par = init.vals[i1, ], 1622 fn = .my.calib.objfunction.rrvglm, 1623 method = optim.control$method.optim, # "BFGS" or "CG" or... 1624 control = c(fnscale = ifelse(minimize.obfunct, 1, -1), 1625 use.optim.control), # as.vector() needed 1626 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 1627 extra = object@extra, 1628 objfun = obfunct, 1629 Object = 777, # object, 1630 Coefs = Coefobject, 1631 B1bix = B1bix[i1, , drop = FALSE], 1632 misc.list = object@misc, 1633 Everything = FALSE, # Differs from Step 5 below 1634 mu.function = mu.function) 1635 1636 if (optim.control$trace) { 1637 if (ans$convergence == 0) 1638 cat("Successful convergence\n") else 1639 cat("Unsuccessful convergence\n") 1640 flush.console() 1641 } 1642 if (ans$convergence == 0) { 1643 BestOFpar[i1, ] <- ans$par 1644 BestOFvalues[i1] <- ans$value 1645 } # else do nothing since NA_real_ signifies convergence failure 1646 } # for i1 1647 1648 1649 1650 1651 1652 1653 1654 prettyCLO <- function(BestOFpar, newdata, Rank) { 1655 if (Rank == 1) { 1656 BestOFpar <- c(BestOFpar) # Drop the dimension 1657 if (!is.null(dimnames(newdata)[[1]])) { 1658 names(BestOFpar) <- dimnames(newdata)[[1]] 1659 } 1660 } else 1661 dimnames(BestOFpar) <- list(dimnames(newdata)[[1]], 1662 param.names("latvar", Rank, skip1 = TRUE)) 1663 BestOFpar 1664 } # prettyCLO 1665 1666 BestOFpar <- prettyCLO(BestOFpar, newdata, Rank) # Dimension may drop. 1667 attr(BestOFpar,"objectiveFunction") <- 1668 prettyCLO(BestOFvalues, newdata, Rank = 1) 1669 if (type == "latvar" && (!cf.confint && !lr.confint)) { 1670 return(BestOFpar) 1671 } 1672 1673 1674 1675 1676 1677 1678 if (lr.confint && Rank > 1) { 1679 warning("argument 'lr.confint' should only be TRUE if Rank == 1. ", 1680 "Setting 'lr.confint = FALSE'.") 1681 lr.confint <- FALSE 1682 } 1683 1684 1685 if (lr.confint && !(type %in% c("latvar", "everything"))) { 1686 warning("argument 'lr.confint' should only be TRUE if ", 1687 "'type = \"latvar\"' or 'type = \"everything\"'. ", 1688 "Setting 'lr.confint = FALSE'.") 1689 lr.confint <- FALSE 1690 } 1691 if (lr.confint && Rank == 1) { 1692 1693 format.perc <- function(probs, digits) 1694 paste(format(100 * probs, trim = TRUE, scientific = FALSE, 1695 digits = digits), "%") 1696 aa <- (1 - level) / 2 1697 aa <- c(aa, 1 - aa) 1698 1699 pct <- format.perc(aa, 3) 1700 cimat <- array(NA, dim = c(nn, 2L), 1701 dimnames = list(dimnames(newdata)[[1]], pct)) 1702 1703 for (i1 in 1:nn) { 1704 if (optim.control$trace) { 1705 cat("Solving for the roots for obsn", i1, "---------------\n") 1706 flush.console() 1707 } 1708 1709 1710 foo3.lhs.rhs <- 1711 function(bnu, 1712 y0, extra = NULL, 1713 objfun, Coefs, 1714 Object, # Not needed 1715 B1bix, 1716 misc.list, 1717 Everything = FALSE, 1718 mu.function, 1719 BestOFvalues = NA, 1720 level = 0.95, 1721 criterion.arg = "loglikelihood") { 1722 if (!(criterion.arg %in% c("deviance", "loglikelihood"))) 1723 stop("'criterion.arg' must be 'deviance' or 'loglikelihood'") 1724 ifelse(criterion.arg == "deviance", 1, -2) * 1725 (-BestOFvalues + 1726 .my.calib.objfunction.rrvglm(bnu = bnu, 1727 y0 = y0, 1728 extra = extra, 1729 objfun = objfun, 1730 Object = Object, 1731 Coefs = Coefs, 1732 B1bix = B1bix, 1733 misc.list = misc.list, 1734 Everything = Everything, 1735 mu.function = mu.function)) - 1736 qchisq(level, df = 1) 1737 } 1738 1739 1740 for (Side in 1:2) { 1741 ans.lhs.rhs <- 1742 uniroot(f = foo3.lhs.rhs, 1743 interval = if (Side == 1) c(Lvec, BestOFpar[i1]) else 1744 c(BestOFpar[i1], Uvec), 1745 extendInt = ifelse(Side == 1, "downX", "upX"), 1746 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 1747 extra = object@extra, 1748 objfun = obfunct, 1749 Object = 777, # object, 1750 Coefs = Coefobject, 1751 B1bix = B1bix[i1, , drop = FALSE], 1752 misc.list = object@misc, 1753 Everything = FALSE, 1754 mu.function = mu.function, 1755 BestOFvalues = BestOFvalues[i1], level = level, 1756 criterion.arg = object@misc$criterion) 1757 cimat[i1, Side] <- ans.lhs.rhs$root 1758 } # Side 1759 } # for i1 1760 1761 if (type == "latvar") 1762 return(cbind(estimate = BestOFpar, 1763 objfun = BestOFvalues, 1764 cimat)) 1765 } # if (lr.confint && Rank == 1) 1766 1767 1768 1769 1770 1771 1772 1773 if (cf.confint && Rank > 1) { 1774 warning("argument 'cf.confint' should only be TRUE if Rank == 1. ", 1775 "Setting 'cf.confint = FALSE'.") 1776 cf.confint <- FALSE 1777 } 1778 1779 1780 if (cf.confint && !(type %in% c("latvar", "everything"))) { 1781 warning("argument 'cf.confint' should only be TRUE if ", 1782 "'type = \"latvar\"' or 'type = \"everything\"'. ", 1783 "Setting 'cf.confint = FALSE'.") 1784 cf.confint <- FALSE 1785 } 1786 if (cf.confint && Rank == 1) { 1787 1788 format.perc <- function(probs, digits) 1789 paste(format(100 * probs, trim = TRUE, scientific = FALSE, 1790 digits = digits), "%") 1791 aa <- (1 - level) / 2 1792 aa <- c(aa, 1 - aa) 1793 pct <- format.perc(aa, 3) 1794 cimat2 <- array(NA, dim = c(nn, 2L), 1795 dimnames = list(dimnames(newdata)[[1]], pct)) 1796 1797 1798 1799 for (i1 in 1:nn) { 1800 if (optim.control$trace) { 1801 cat("\nSolving for the roots for obsn", i1, "---------------\n") 1802 flush.console() 1803 } 1804 1805 1806 foo4.lhs.rhs <- 1807 function(bnu, 1808 y0, extra = NULL, 1809 objfun, Coefs, 1810 Object, # Not needed 1811 B1bix, 1812 misc.list, 1813 Everything = FALSE, 1814 mu.function, 1815 BestOFvalues = NA, 1816 pr.level = c(0.05, 0.95)[1] 1817 ) { 1818 charfun.clo.cdf(bnu = bnu, 1819 y0 = y0, 1820 extra = extra, 1821 objfun = objfun, 1822 Object = Object, 1823 Coefs = Coefs, 1824 B1bix = B1bix, 1825 misc.list = misc.list, 1826 Everything = Everything, 1827 mu.function = mu.function 1828 ) - 1829 pr.level 1830 } 1831 1832 1833 1834 1835 1836 1837 1838 for (Side in 1:2) { 1839 ans.lhs.rhs <- 1840 uniroot(f = foo4.lhs.rhs, 1841 interval = if (Side == 1) c(Lvec, BestOFpar[i1]) else 1842 c(BestOFpar[i1], Uvec), 1843 extendInt = "yes", # Might be worse than above. 1844 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 1845 extra = object@extra, 1846 objfun = obfunct, 1847 Object = object, 1848 Coefs = Coefobject, 1849 B1bix = B1bix[i1, , drop = FALSE], 1850 misc.list = object@misc, 1851 Everything = FALSE, 1852 mu.function = mu.function, 1853 BestOFvalues = BestOFvalues[i1], 1854 pr.level = ifelse(Side == 1, aa[1], aa[2]) 1855 ) 1856 cimat2[i1, Side] <- ans.lhs.rhs$root 1857 } # Side 1858 } # for i1 1859 1860 vecTF <- cimat2[, 2] < cimat2[, 1] 1861 if (any(vecTF)) { 1862 temp <- cimat2[vecTF, 1] 1863 cimat2[vecTF, 1] <- cimat2[vecTF, 2] 1864 cimat2[vecTF, 2] <- temp 1865 } 1866 if (type == "latvar") 1867 return(cbind(estimate = BestOFpar, 1868 objfun = BestOFvalues, 1869 cimat2)) 1870 1871 } # if (cf.confint && Rank == 1) 1872 1873 1874 1875 1876 1877 1878 1879 1880 etaValues <- matrix(NA_real_, nn, M) 1881 muValues <- matrix(NA_real_, nn, ncol(fitted(object))) 1882 vcValues <- if (get.SEs) array(0, c(Rank, Rank, nn)) else NULL 1883 1884 1885 if (optim.control$trace) 1886 cat("\n") 1887 1888 for (i1 in 1:nn) { 1889 if (optim.control$trace) { 1890 cat("Evaluating quantities for observation", i1, 1891 "-----------------\n") 1892 flush.console() 1893 } 1894 ans5 <- .my.calib.objfunction.rrvglm( 1895 bnu = if (Rank == 1) BestOFpar[i1] else BestOFpar[i1, ], 1896 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 1897 extra = object@extra, 1898 objfun = obfunct, # deviance 1899 Object = 777, # object, 1900 Coefs = Coefobject, 1901 B1bix = B1bix[i1, , drop = FALSE], 1902 misc.list = object@misc, 1903 Everything = TRUE, # Differs from Step 3. 1904 mu.function = mu.function) 1905 1906 1907 muValues[i1, ] <- ans5$mu 1908 etaValues[i1, ] <- ans5$eta 1909 1910 1911 if (get.SEs) 1912 vcValues[, , i1] <- if (se.type == "dzwald") 1913 dzwald.rrvglm( 1914 bnu0 = if (Rank == 1) BestOFpar[i1] else BestOFpar[i1, ], 1915 y0 = newdata[i1, , drop = FALSE], # drop added 20150624 1916 Object = object, 1917 CoefsObject = Coefobject, 1918 B1bix = B1bix[i1, , drop = FALSE], 1919 mu.function = mu.function 1920 ) else 1921 ans5$vcmat 1922 1923 1924 1925 } # for i1 1926 1927 1928 1929 1930 1931 1932 1933 dimnames(muValues) <- dimnames(newdata) 1934 dimnames(etaValues) <- list(dimnames(newdata)[[1]], 1935 dimnames(object@predictors)[[2]]) 1936 if (get.SEs) 1937 dimnames(vcValues) <- list(as.character(1:Rank), 1938 as.character(1:Rank), 1939 dimnames(newdata)[[1]]) 1940 1941 switch(type, 1942 latvar = BestOFpar, # Done already, so not really needed 1943 predictors = etaValues, 1944 response = muValues, 1945 vcov = vcValues, 1946 everything = list( 1947 latvar = BestOFpar, 1948 predictors = etaValues, 1949 response = muValues, 1950 vcov = vcValues, 1951 lr.confint = if (lr.confint) 1952 cbind(estimate = BestOFpar, 1953 objfun = BestOFvalues, 1954 cimat) else NULL, 1955 cf.confint = if (cf.confint) 1956 cbind(estimate = BestOFpar, 1957 objfun = BestOFvalues, 1958 cimat2) else NULL) 1959 ) 1960} # calibrate.rrvglm 1961 1962 1963 1964 1965 1966 1967.my.calib.objfunction.rrvglm <- 1968 function(bnu, 1969 y0, extra = NULL, 1970 objfun, Coefs, 1971 Object, # Not needed 1972 B1bix, 1973 misc.list, 1974 Everything = TRUE, 1975 mu.function) { 1976 1977 1978 1979 bnumat <- cbind(bnu) 1980 Rank <- length(bnu) 1981 eta <- cbind(c(B1bix)) + Coefs@A %*% bnumat 1982 M <- misc.list$M 1983 eta <- matrix(eta, 1, M, byrow = TRUE) 1984 1985 mu <- matrix(mu.function(eta, extra = extra), nrow = 1) 1986 obvalue <- objfun(mu = mu, y = y0, 1987 w = 1, # ignore prior.weights on the object zz 1988 residuals = FALSE, eta = eta, extra = extra) 1989 if (Everything) { 1990 vcmat <- diag(Rank) 1991 vcmat <- solve(vcmat) 1992 } else { 1993 vcmat <- NULL 1994 } 1995 if (Everything) 1996 list(eta = eta, 1997 mu = mu, 1998 obvalue = obvalue, 1999 vcmat = vcmat) else 2000 obvalue 2001} # .my.calib.objfunction.rrvglm 2002 2003 2004 2005 2006 2007 2008calibrate.rrvglm.control <- 2009 function(object, 2010 trace = FALSE, # passed into optim() 2011 method.optim = "BFGS", # passed into optim(method = Method) 2012 gridSize = ifelse(Rank == 1, 17, 9), 2013 ...) { 2014 2015 2016 Rank <- object@control$Rank 2017 if (!is.Numeric(gridSize, positive = TRUE, 2018 integer.valued = TRUE, length.arg = 1)) 2019 stop("bad input for argument 'gridSize'") 2020 if (gridSize < 2) 2021 stop("argument 'gridSize' must be >= 2") 2022 2023 list( 2024 trace = as.numeric(trace)[1], 2025 method.optim = method.optim, 2026 gridSize = gridSize 2027 ) 2028} # calibrate.rrvglm.control 2029 2030 2031 2032 2033 2034setMethod("calibrate", "rrvglm", function(object, ...) 2035 calibrate.rrvglm(object, ...)) 2036 2037 2038 2039 setMethod("calibrate", "qrrvglm", function(object, ...) 2040 calibrate.qrrvglm(object, ...)) 2041 2042 2043 2044