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