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 rcim <-
20  function(y,
21           family = poissonff,
22           Rank = 0,
23           M1 = NULL,
24           weights = NULL,
25           which.linpred = 1,
26           Index.corner = ifelse(is.null(str0), 0, max(str0)) + 1:Rank,
27           rprefix = "Row.",
28           cprefix = "Col.",
29           iprefix = "X2.",
30           offset = 0,
31
32           str0 = if (Rank) 1 else NULL,  # Ignored if Rank == 0
33           summary.arg = FALSE, h.step = 0.0001,
34           rbaseline = 1, cbaseline = 1,
35
36           has.intercept = TRUE,
37
38           M = NULL,
39
40           rindex = 2:nrow(y),  # Row index
41           cindex = 2:ncol(y),  # Col index
42           iindex = 2:nrow(y),  # Interaction index
43
44
45
46           ...) {
47
48
49
50
51
52  rindex <- unique(sort(rindex))
53  cindex <- unique(sort(cindex))
54  iindex <- unique(sort(iindex))
55
56
57  if (Rank == 0 && !has.intercept)
58    warning("probably 'has.intercept == TRUE' is better for ",
59            "a rank-0 model")
60
61
62
63  ncoly <- ncol(y)
64
65
66  noroweffects <- FALSE
67  nocoleffects <- FALSE
68
69  if (!is.Numeric(which.linpred, length.arg = 1,
70                  integer.valued = TRUE, positive = TRUE))
71    stop("bad input for argument 'which.linpred'")
72
73  if (!is.character(rprefix))
74    stop("argument 'rprefix' must be character")
75  if (!is.character(cprefix))
76    stop("argument 'cprefix' must be character")
77
78  if (is.character(family))
79    family <- get(family)
80  if (is.function(family))
81    family <- ((family)())
82  if (!inherits(family, "vglmff")) {
83    stop("'family = ", family, "' is not a VGAM family function")
84  }
85  efamily <- family
86
87
88  if (!is.Numeric(M1)) {
89    iefamily <- efamily@infos
90
91    if (is.function(iefamily))
92      M1 <- (iefamily())$M1
93      if (is.Numeric(M1))
94        M1 <- abs(M1)
95  }
96  if (!is.Numeric(M1)) {
97    if (!is.Numeric(M))
98      warning("cannot determine the value of 'M1'.",
99              "Assuming the value one.")
100    M1 <- 1
101  }
102
103
104
105  M <- if (is.null(M)) M1 * ncol(y) else M
106
107  special <- (M > 1) && (M1 == 1)
108
109
110
111  object.save <- y
112  y <- if (is(y, "rrvglm")) {
113    depvar(object.save)
114  } else {
115    as(as.matrix(y), "matrix")
116  }
117  if (length(dim(y)) != 2 || nrow(y) < 3 || ncol(y) < 3)
118    stop("argument 'y' must be a matrix with >= 3 rows & columns, or ",
119         "a rrvglm() object")
120
121
122  .rcim.df <-
123    if (!noroweffects)
124      data.frame("Row.2" = I.col(2, nrow(y))) else  # See below
125    if (!nocoleffects)
126      data.frame("Col.2" = I.col(2, nrow(y))) else  # See below
127      stop("at least one of 'noroweffects' and 'nocoleffects' ",
128           "must be FALSE")
129
130
131  min.row.val <- rindex[1]  # == min(rindex) since its sorted # Usually 2
132  min.col.val <- cindex[1]  # == min(cindex) since its sorted # Usually 2
133  if (!noroweffects) {
134    colnames( .rcim.df ) <-
135      paste(rprefix, as.character(min.row.val),  # "2",
136            sep = "")  # Overwrite "Row.2"
137  } else if (!nocoleffects) {
138    colnames( .rcim.df ) <-
139      paste(cprefix, as.character(min.col.val),  # "2",
140            sep = "")  # Overwrite "Col.2"
141  }
142
143
144
145  yn1 <- if (length(dimnames(y)[[1]])) dimnames(y)[[1]] else
146             param.names(iprefix, nrow(y))
147  warn.save <- options()$warn
148  options(warn = -3)  # Suppress the warnings (hopefully, temporarily)
149  if (any(!is.na(as.numeric(substring(yn1, 1, 1)))))
150    yn1 <- param.names(iprefix, nrow(y))
151  options(warn = warn.save)
152
153
154  nrprefix <- as.name(rprefix)
155  ncprefix <- as.name(cprefix)
156
157
158  assign(rprefix, factor(1:nrow(y)))
159  modmat.row <- substitute(
160           model.matrix( ~ .rprefix ), list( .rprefix = nrprefix ))
161
162  LLL <- ifelse(special, M, ncol(y))
163  assign(cprefix, factor(1:LLL))
164
165
166  modmat.col <- substitute(
167           model.matrix( ~ .cprefix ), list( .cprefix = ncprefix ))
168  modmat.row <- eval( modmat.row )
169  modmat.col <- eval( modmat.col )
170
171
172
173
174  Hlist <-
175    if (has.intercept) {
176      list("(Intercept)" = matrix(1, LLL, 1))
177    } else {
178      temp <- list("Row.2" = matrix(1, LLL, 1))  # Overwrite this name:
179      names(temp) <- paste(rprefix, as.character(min.row.val), sep = "")
180      temp
181    }
182
183
184  if (!noroweffects)
185    for (ii in rindex) {
186         Hlist[[paste(rprefix, ii, sep = "")]] <- matrix(1, LLL, 1)
187      .rcim.df[[paste(rprefix, ii, sep = "")]] <- modmat.row[, ii]
188    }
189
190
191  if (!nocoleffects)
192    for (ii in cindex) {
193      temp6.mat <- modmat.col[, ii, drop = FALSE]
194         Hlist[[paste(cprefix, ii, sep = "")]] <- temp6.mat
195      .rcim.df[[paste(cprefix, ii, sep = "")]] <- rep_len(1, nrow(y))
196    }
197
198
199  if (Rank > 0) {
200    for (ii in iindex) {
201
202      Hlist[[yn1[ii]]] <- diag(LLL)
203
204      .rcim.df[[yn1[ii]]] <- I.col(ii, nrow(y))
205    }
206  }
207
208
209  dimnames(.rcim.df) <- list(if (length(dimnames(y)[[1]]))
210                               dimnames(y)[[1]] else
211                               as.character(iindex),
212                             dimnames( .rcim.df )[[2]])
213
214  str1 <- paste(if (has.intercept) "~ 1 + " else "~ -1 + ", rprefix,
215                as.character(min.row.val),  # "2",
216                sep = "")
217
218
219  if (nrow(y) > 2)
220    str1 <- paste(str1,
221                  paste(rprefix, rindex[-1], sep = "", collapse = " + "),
222                  sep = " + ")
223
224
225    str1 <- paste(str1,
226                  paste(cprefix, cindex, sep = "", collapse = " + "),
227                  sep = " + ")
228
229
230
231
232  str2 <- paste("y", str1)
233  if (Rank > 0) {
234    str2 <- paste(str2,
235                  paste(yn1[iindex], sep = "", collapse = " + "),
236                  sep = " + ")
237  }
238
239
240
241
242  controlfun <- if (Rank == 0)   vglm.control else rrvglm.control  # orig.
243
244
245  mycontrol <- controlfun(Rank = Rank,
246                          Index.corner = Index.corner,
247                          str0 = str0, ...)
248
249  if (mycontrol$trace) {
250  }
251
252
253
254  if ((mindim <- min(nrow(y), ncol(y))) <= Rank) {
255    stop("argument 'Rank' is too high. Must be a value from 0 ",
256         "to ", mindim - 1, " inclusive")
257  }
258
259
260
261  if (Rank > 0)
262    mycontrol$noRRR <- as.formula(str1)  # Overwrite this
263
264  assign(".rcim.df", .rcim.df, envir = VGAM::VGAMenv)
265
266  warn.save <- options()$warn
267  options(warn = -3)  # Suppress the warnings (hopefully, temporarily)
268
269  if (mycontrol$trace) {
270  }
271
272
273  if (M1 > 1) {
274    orig.Hlist <- Hlist
275
276    kmat1 <- matrix(0, nrow = M1, ncol = 1)
277    kmat1[which.linpred, 1] <- 1
278    kmat0 <- (diag(M1))[, -which.linpred, drop = FALSE]
279
280    for (ii in seq_along(Hlist)) {
281      Hlist[[ii]] <- kronecker(Hlist[[ii]], kmat1)
282    }
283    if (has.intercept)
284      Hlist[["(Intercept)"]] <- cbind(Hlist[["(Intercept)"]],
285                                      kronecker(matrix(1, ncoly, 1),
286                                                kmat0))
287
288
289    if (mycontrol$trace) {
290    }
291  }
292
293
294
295  offset.matrix <-
296    matrix(offset, nrow = nrow(y),
297                   ncol = M)  # byrow = TRUE
298
299
300
301  answer <- if (Rank > 0) {
302    if (is(object.save, "rrvglm")) object.save else
303      rrvglm(as.formula(str2),
304             family = family,
305             constraints = Hlist,
306             offset = offset.matrix,
307             weights = if (length(weights))
308                       weights else rep_len(1, nrow(y)),
309             ...,
310             control = mycontrol, data = .rcim.df )
311  } else {
312    if (is(object.save, "vglm")) object.save else
313        vglm(as.formula(str2),
314             family = family,
315             constraints = Hlist,
316             offset = offset.matrix,
317             weights = if (length(weights))
318                       weights else rep_len(1, nrow(y)),
319             ...,
320             control = mycontrol, data = .rcim.df )
321  }
322
323  options(warn = warn.save)  # Restore warnings back to prior state
324
325
326  answer <- if (summary.arg) {
327    if (Rank > 0) {
328      summary.rrvglm(as(answer, "rrvglm"), h.step = h.step)
329    } else {
330      summary(answer)
331    }
332  } else {
333    as(answer, ifelse(Rank > 0, "rcim",  "rcim0"))
334  }
335
336
337  answer@misc$rbaseline     <- rbaseline
338  answer@misc$cbaseline     <- cbaseline
339  answer@misc$which.linpred <- which.linpred
340  answer@misc$offset        <- offset.matrix
341
342  answer
343}
344
345
346
347
348
349
350
351
352summaryrcim <- function(object, ...) {
353  rcim(depvar(object), summary.arg = TRUE, ...)
354}
355
356
357
358
359
360
361
362
363
364 setClass("rcim0", representation(not.needed = "numeric"),
365          contains = "vglm")  # Added 20110506
366
367 setClass("rcim", representation(not.needed = "numeric"),
368          contains = "rrvglm")
369
370
371setMethod("summary", "rcim0",
372          function(object, ...)
373          summaryrcim(object, ...))
374
375
376setMethod("summary", "rcim",
377          function(object, ...)
378          summaryrcim(object, ...))
379
380
381
382
383
384
385
386
387
388
389 Rcim <- function(mat, rbaseline = 1, cbaseline = 1) {
390
391  mat <- as.matrix(mat)
392  RRR <- dim(mat)[1]
393  CCC <- dim(mat)[2]
394
395  rnames <- if (is.null(rownames(mat))) {
396    param.names("X", RRR)
397  } else {
398    rownames(mat)
399  }
400
401  cnames <- if (is.null(colnames(mat))) {
402    param.names("Y", CCC)
403  } else {
404    colnames(mat)
405  }
406
407  r.index <- if (is.character(rbaseline))
408               which(rownames(mat) == rbaseline) else
409                     if (is.numeric(rbaseline)) rbaseline else
410                         stop("argement 'rbaseline' must be numeric",
411                               "or character of the level of row")
412
413  c.index <- if (is.character(cbaseline))
414               which(colnames(mat) == cbaseline) else
415                     if (is.numeric(cbaseline)) cbaseline else
416                         stop("argement 'cbaseline' must be numeric",
417                               "or character of the level of row")
418
419  if (length(r.index) != 1)
420    stop("Could not match with argument 'rbaseline'")
421
422  if (length(c.index) != 1)
423    stop("Could not match with argument 'cbaseline'")
424
425
426  yswap <- rbind(mat[r.index:RRR, ],
427                 if (r.index > 1) mat[1:(r.index - 1),] else NULL)
428  yswap <- cbind(yswap[, c.index:CCC],
429                 if (c.index > 1) yswap[, 1:(c.index - 1)] else NULL)
430
431  new.rnames <- rnames[c(r.index:RRR,
432                         if (r.index > 1) 1:(r.index - 1) else NULL)]
433  new.cnames <- cnames[c(c.index:CCC,
434                         if (c.index > 1) 1:(c.index - 1) else NULL)]
435  colnames(yswap) <- new.cnames
436  rownames(yswap) <- new.rnames
437
438  yswap
439}
440
441
442
443
444
445
446
447
448
449
450
451
452
453 plotrcim0  <- function(object,
454     centered = TRUE, which.plots = c(1, 2),
455     hline0 = TRUE, hlty = "dashed", hcol = par()$col, hlwd = par()$lwd,
456     rfirst = 1, cfirst = 1,
457     rtype = "h", ctype = "h",
458     rcex.lab = 1, rcex.axis = 1,  # rlabels = FALSE,
459     rtick = FALSE,
460     ccex.lab = 1, ccex.axis = 1,  # clabels = FALSE,
461     ctick = FALSE,
462     rmain = "Row effects", rsub = "",
463     rxlab = "", rylab = "Row effects",
464     cmain = "Column effects", csub = "",
465     cxlab = "", cylab = "Column effects",
466     rcol = par()$col, ccol = par()$col,
467     no.warning = FALSE,
468     ...) {
469
470
471  nparff <- if (is.numeric(object@family@infos()$M1)) {
472    object@family@infos()$M1
473  } else {
474    1
475  }
476
477
478
479  if (!no.warning &&
480      is.numeric(object@control$Rank) &&
481      object@control$Rank != 0)
482    warning("argument 'object' is not Rank-0")
483
484
485  n.lm  <- nrow(object@y)
486
487  cobj <- coefficients(object)
488
489  upperbound <- if (!is.numeric(object@control$Rank) ||
490                   object@control$Rank == 0) length(cobj) else
491               length(object@control$colx1.index)
492
493  orig.roweff <- c("Row.1" = 0, cobj[(nparff + 1) : (nparff + n.lm - 1)])
494  orig.coleff <- c("Col.1" = 0, cobj[(nparff + n.lm) : upperbound])
495  last.r <- length(orig.roweff)
496  last.c <- length(orig.coleff)
497
498
499  orig.raxisl  <- rownames(object@y)
500  orig.caxisl  <- colnames(object@y)
501  if (is.null(orig.raxisl))
502    orig.raxisl <- as.character(1:nrow(object@y))
503  if (is.null(orig.caxisl))
504    orig.caxisl <- as.character(1:ncol(object@y))
505
506  roweff.orig <-
507  roweff <- orig.roweff[c(rfirst:last.r,
508                          if (rfirst > 1) 1:(rfirst-1) else NULL)]
509  coleff.orig <-
510  coleff <- orig.coleff[c(cfirst:last.c,
511                          if (cfirst > 1) 1:(cfirst-1) else NULL)]
512
513  if (centered) {
514    roweff <- scale(roweff, scale = FALSE)  # Center it only
515    coleff <- scale(coleff, scale = FALSE)  # Center it only
516  }
517
518  raxisl <- orig.raxisl[c(rfirst:last.r,
519                          if (rfirst > 1) 1:(rfirst-1) else NULL)]
520
521  caxisl <- orig.caxisl[c(cfirst:last.c,
522                          if (cfirst > 1) 1:(cfirst-1) else NULL)]
523
524
525  if (any(which.plots == 1, na.rm = TRUE)) {
526    plot(roweff, type = rtype,
527         axes = FALSE, col = rcol, main = rmain,
528         sub  = rsub, xlab = rxlab, ylab = rylab, ...)
529
530    axis(1, at = seq_along(raxisl),
531         cex.lab = rcex.lab,
532         cex.axis = rcex.axis,
533         labels = raxisl)
534    axis(2, cex.lab = rcex.lab, ...)  # las = rlas)
535
536    if (hline0)
537      abline(h = 0, lty = hlty, col = hcol, lwd = hlwd)
538  }
539
540
541  if (any(which.plots == 2, na.rm = TRUE)) {
542    plot(coleff, type = ctype,
543         axes = FALSE, col = ccol, main = cmain,  # lwd = 2, xpd = FALSE,
544         sub  = csub, xlab = cxlab, ylab = cylab, ...)
545
546    axis(1, at = seq_along(caxisl),
547         cex.lab = ccex.lab,
548         cex.axis = ccex.axis,
549         labels = caxisl)
550    axis(2, cex.lab = ccex.lab, ...)  # las = clas)
551
552    if (hline0)
553      abline(h = 0, lty = hlty, col = hcol, lwd = hlwd)
554  }
555
556
557
558
559  object@post$row.effects = roweff
560  object@post$col.effects = coleff
561  object@post$raw.row.effects = roweff.orig
562  object@post$raw.col.effects = coleff.orig
563
564  invisible(object)
565}
566
567
568
569
570
571setMethod("plot", "rcim0",
572          function(x, y, ...)
573          plotrcim0(object = x, ...))
574
575
576setMethod("plot", "rcim",
577          function(x, y, ...)
578          plotrcim0(object = x, ...))
579
580
581
582
583
584
585
586
587
588
589
590
591
592 moffset <-
593  function(mat, roffset = 0, coffset = 0, postfix = "",
594           rprefix = "Row.",
595           cprefix = "Col."
596          ) {
597
598
599
600
601
602  if ((is.numeric(roffset) && (roffset == 0)) &&
603      (is.numeric(coffset) && (coffset == 0)))
604    return(mat)
605
606
607  vecmat <- c(unlist(mat))
608  ind1 <- if (is.character(roffset))
609             which(rownames(mat) == roffset) else
610                   if (is.numeric(roffset)) roffset + 1 else
611                     stop("argument 'roffset' not matched (character).",
612                           " It must be numeric, ",
613                           "else character and match the ",
614                           "row names of the response")
615  ind2 <- if (is.character(coffset))
616             which(colnames(mat) == coffset) else
617                   if (is.numeric(coffset)) coffset + 1 else
618                     stop("argument 'coffset' not matched (character).",
619                           " It must be numeric, ",
620                           "else character and match the ",
621                           "column names of the response")
622
623  if (!is.Numeric(ind1, positive = TRUE,
624                  integer.valued = TRUE, length.arg = 1) ||
625      !is.Numeric(ind2, positive = TRUE,
626                  integer.valued = TRUE, length.arg = 1))
627    stop("bad input for arguments 'roffset' and/or 'coffset'")
628  if (ind1 > nrow(mat))
629    stop("too large a value for argument 'roffset'")
630  if (ind2 > ncol(mat))
631    stop("too large a value for argument 'coffset'")
632
633
634  start.ind <- (ind2 - 1)* nrow(mat) + ind1
635
636
637  svecmat <- vecmat[c(start.ind:(nrow(mat) * ncol(mat)),
638                     0:(start.ind - 1))]
639
640  rownames.mat <- rownames(mat)
641  if (length(rownames.mat) != nrow(mat))
642    rownames.mat <- param.names(rprefix, nrow(mat))
643
644  colnames.mat <- colnames(mat)
645  if (length(colnames.mat) != ncol(mat))
646    colnames.mat <- param.names(cprefix, ncol(mat))
647
648
649  newrn <- if (roffset > 0)
650            c(rownames.mat[c(ind1:nrow(mat))],
651              paste(rownames.mat[0:(ind1-1)], postfix, sep = "")) else
652           rownames.mat
653
654  newcn <- c(colnames.mat[c(ind2:ncol(mat), 0:(ind2 - 1))])
655  if (roffset > 0)
656    newcn <- paste(newcn, postfix, sep = "")
657
658  newmat <- matrix(svecmat, nrow(mat), ncol(mat),
659                  dimnames = list(newrn, newcn))
660  newmat
661}
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681Confint.rrnb <- function(rrnb2, level = 0.95) {
682
683  if (class(rrnb2) != "rrvglm")
684    stop("argument 'rrnb2' does not appear to be a rrvglm() object")
685
686  if (!any(rrnb2@family@vfamily == "negbinomial"))
687    stop("argument 'rrnb2' does not appear to be a negbinomial() fit")
688
689  if (rrnb2@control$Rank != 1)
690    stop("argument 'rrnb2' is not Rank-1")
691
692  if (rrnb2@misc$M != 2)
693    stop("argument 'rrnb2' does not have M = 2")
694
695  if (!all(rrnb2@misc$link == "loglink"))
696    stop("argument 'rrnb2' does not have log links for both parameters")
697
698  a21.hat <- (Coef(rrnb2)@A)["loglink(size)", 1]
699  beta11.hat <- Coef(rrnb2)@B1["(Intercept)", "loglink(mu)"]
700  beta21.hat <- Coef(rrnb2)@B1["(Intercept)", "loglink(size)"]
701  delta1.hat <- exp(a21.hat * beta11.hat - beta21.hat)
702  delta2.hat <- 2 - a21.hat
703
704  SE.a21.hat <- sqrt(vcovrrvglm(rrnb2)["I(latvar.mat)", "I(latvar.mat)"])
705
706
707  ci.a21 <- a21.hat +  c(-1, 1) * qnorm(1 - (1 - level)/2) * SE.a21.hat
708  (ci.delta2 <- 2 - rev(ci.a21))  # e.g., the 95 percent CI
709
710  list(a21.hat    = a21.hat,
711       beta11.hat = beta11.hat,
712       beta21.hat = beta21.hat,
713       CI.a21     = ci.a21,
714       CI.delta2  = ci.delta2,
715       delta1     = delta1.hat,
716       delta2     = delta2.hat,
717       SE.a21.hat = SE.a21.hat)
718}
719
720
721
722
723
724Confint.nb1 <- function(nb1, level = 0.95) {
725
726
727
728
729  if (class(nb1) != "vglm")
730    stop("argument 'nb1' does not appear to be a vglm() object")
731
732  if (!any(nb1@family@vfamily == "negbinomial"))
733    stop("argument 'nb1' does not appear to be a negbinomial() fit")
734
735  if (!all(unlist(constraints(nb1)[-1]) == 1))
736    stop("argument 'nb1' does not appear to have 'parallel = TRUE'")
737
738  if (!all(unlist(constraints(nb1)[1]) == c(diag(nb1@misc$M))))
739    stop("argument 'nb1' does not have 'parallel = FALSE' ",
740         "for the intercept")
741
742  if (nb1@misc$M != 2)
743    stop("argument 'nb1' does not have M = 2")
744
745  if (!all(nb1@misc$link == "loglink"))
746    stop("argument 'nb1' does not have log links for both parameters")
747
748  cnb1 <- coefficients(as(nb1, "vglm"), matrix = TRUE)
749  mydiff <- (cnb1["(Intercept)", "loglink(size)"] -
750             cnb1["(Intercept)", "loglink(mu)"])
751  delta0.hat <- exp(mydiff)
752  (phi0.hat <- 1 + 1 / delta0.hat)  # MLE of phi0
753
754
755  myvcov <- vcov(as(nb1, "vglm"))  # Not great; improve this!
756
757
758
759  myvec <- cbind(c(-1, 1, rep_len(0, nrow(myvcov) - 2)))
760  se.mydiff <- c(sqrt(t(myvec) %*% myvcov %*% myvec))  # c() to drop the 2-D
761
762  ci.mydiff <- mydiff + c(-1, 1) * qnorm(1 - (1 - level)/2) * se.mydiff
763
764  ci.delta0 <- ci.exp.mydiff <- exp(ci.mydiff)
765  (ci.phi0 <- 1 + 1 / rev(ci.delta0))  # e.g., the 95 percent CI for phi0
766
767  list(CI.phi0    = ci.phi0,
768       CI.delta0  = ci.delta0,
769       delta0     = delta0.hat,
770       phi0       = phi0.hat)
771}
772
773
774
775
776
777
778plota21 <- function(rrvglm2, show.plot = TRUE, nseq.a21 = 31,
779                    se.eachway = c(5, 5),  # == c(LHS, RHS),
780                    trace.arg = TRUE,
781                    lwd = 2, ...) {
782
783
784
785
786
787  if (class(rrvglm2) != "rrvglm")
788    stop("argument 'rrvglm2' does not appear to be a rrvglm() object")
789
790  if (rrvglm2@control$Rank != 1)
791    stop("argument 'rrvglm2' is not Rank-1")
792
793  if (rrvglm2@misc$M != 2)
794    stop("argument 'rrvglm2' does not have M = 2")
795
796
797  loglik.orig <- logLik(rrvglm2)
798  temp1 <- Confint.rrnb(rrvglm2)  # zz
799
800  a21.hat <- (Coef(rrvglm2)@A)[2, 1]
801  SE.a21.hat <- temp1$SE.a21.hat
802
803
804  SE.a21.hat <- sqrt(vcov(rrvglm2)["I(latvar.mat)", "I(latvar.mat)"])
805
806
807  big.ci.a21 <- a21.hat +  c(-1, 1) * se.eachway * SE.a21.hat
808  seq.a21 <- seq(big.ci.a21[1], big.ci.a21[2], length = nseq.a21)
809  Hlist.orig <- constraints.vlm(rrvglm2, type = "term")
810
811
812  alreadyComputed <- !is.null(rrvglm2@post$a21.matrix)
813
814
815  a21.matrix <- if (alreadyComputed) rrvglm2@post$a21.matrix else
816                cbind(a21 = seq.a21, loglikelihood = 0)
817  prev.etastart <- predict(rrvglm2)  # Halves the computing time
818  funname <- "vglm"
819  listcall <- as.list(rrvglm2@call)
820
821
822  if (!alreadyComputed)
823  for (ii in 1:nseq.a21) {
824    if (trace.arg)
825      print(ii)
826    argslist <- vector("list", length(listcall) - 1)
827    for (kay in 2:(length(listcall)))
828      argslist[[kay - 1]] <- listcall[[kay]]
829
830    names(argslist) <- c(names(listcall)[-1])
831
832    argslist$trace       <- trace.arg
833    argslist$etastart    <- prev.etastart
834    argslist$constraints <- Hlist.orig
835
836
837    for (kay in 2:length(argslist[["constraints"]])) {
838       argslist[["constraints"]][[kay]] <- rbind(1, a21.matrix[ii, 1])
839    }
840
841
842    fitnew <- do.call(what = funname, args = argslist)
843
844    a21.matrix[ii, 2] <- logLik(fitnew)
845
846    prev.etastart <- predict(fitnew)
847  }
848
849
850
851  if (show.plot) {
852    plot(a21.matrix[ ,1], a21.matrix[ ,2], type = "l",
853            col = "blue", cex.lab = 1.1,
854            xlab = expression(a[21]), ylab = "Log-likelihood")  # ...
855
856    abline(v = (Hlist.orig[[length(Hlist.orig)]])[2, 1],
857           col = "darkorange", lty = "dashed")
858
859    abline(h = loglik.orig,
860           col = "darkorange", lty = "dashed")
861
862    abline(h = loglik.orig -
863               qchisq(0.95, df = 1) / 2,
864           col = "darkorange", lty = "dashed")
865
866    abline(v = a21.hat +  c(-1, 1) * 1.96 * SE.a21.hat,
867           col = "gray50", lty = "dashed", lwd = lwd)
868
869  }  # End of (show.plot)
870
871  rrvglm2@post <- list(a21.matrix = a21.matrix)
872  invisible(rrvglm2)
873}
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888 Qvar <- function(object,
889                  factorname = NULL,
890                  which.linpred = 1,
891                  coef.indices = NULL,
892                  labels = NULL, dispersion = NULL,
893                  reference.name = "(reference)",
894                  estimates = NULL
895                 ) {
896
897
898
899
900
901
902
903
904  if (!is.Numeric(which.linpred, length.arg = 1,
905                  integer.valued = TRUE, positive = TRUE))
906    stop("argument 'which.linpred' must be a positive integer")
907
908
909
910  coef.indices.saved <- coef.indices
911
912  if (!is.matrix(object)) {
913    model <- object
914    if (is.null(factorname) && is.null(coef.indices)) {
915      stop("arguments 'factorname' and 'coef.indices' are ",
916           "both NULL")
917    }
918
919
920    if (is.null(coef.indices)) {
921
922      M <- npred(model)
923      if (M < which.linpred)
924        stop("argument 'which.linpred' must be a value from the set 1:",
925             M)
926
927
928      newfactorname <- if (M > 1) {
929        clist <- constraints(model, type = "term")
930
931        Hk <- clist[[factorname]]
932        Mdot <- ncol(Hk)
933        Hk.row <- Hk[which.linpred, ]
934        if (sum(Hk.row != 0) > 1)
935          stop("cannot handle rows of constraint matrices with more ",
936               "than one nonzero value")
937
938        foo <- function(ii)
939          switch(as.character(ii), '1'="1st", '2'="2nd", '3'="3rd",
940                 paste(ii, "th", sep = ""))
941        if (sum(Hk.row != 0) == 0)
942          stop("factor '", factorname, "' is not used the ",
943               foo(which.linpred), " eta (linear predictor)")
944
945        row.index <- (1:Mdot)[Hk.row != 0]
946
947        all.labels <- vlabel(factorname, ncolHlist = Mdot, M = M)
948        all.labels[row.index]
949      } else {
950        factorname
951      }
952
953      colptr <- attr(model.matrix(object, type = "vlm"), "vassign")
954      colptr <- if (M > 1) {
955        colptr[newfactorname]
956      } else {
957        colptr[[newfactorname]]
958      }
959
960      coef.indices <- colptr
961
962      contmat <- if (length(model@xlevels[[factorname]]) ==
963                     length(coef.indices)) {
964        diag(length(coef.indices))
965      } else {
966        eval(call(model@contrasts[[factorname]],
967                  model@xlevels  [[factorname]]))
968      }
969      rownames(contmat) <- model@xlevels[[factorname]]
970
971      if (is.null(estimates)) {
972        if (M > 1) {
973          estimates <- matrix(-1, nrow(contmat), 1)  # Used to be nc = Mdot
974          ii <- 1
975          estimates[, ii] <- contmat %*%
976                             (coefvlm(model)[(coef.indices[[ii]])])
977        } else {
978          estimates <- contmat %*% (coefvlm(model)[coef.indices])
979        }
980      }
981
982
983
984
985      Covmat <- vcov(model, dispersion = dispersion)
986
987
988
989      covmat <- Covmat[unlist(coef.indices),
990                       unlist(coef.indices), drop = FALSE]
991      covmat <- if (M > 1) {
992        ii <- 1
993        ans <- contmat %*% Covmat[(colptr[[ii]]),
994                                  (colptr[[ii]])] %*% t(contmat)
995        ans
996      } else {
997        contmat %*% covmat %*% t(contmat)
998      }
999    } else { # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
1000      kk <- length(coef.indices)
1001      refPos <- numeric(0)
1002      if (0 %in% coef.indices) {
1003        refPos <- which(coef.indices == 0)
1004        coef.indices <- coef.indices[-refPos]
1005      }
1006
1007
1008
1009
1010      covmat <- vcov(model, dispersion = dispersion)
1011
1012
1013
1014
1015      covmat <- covmat[coef.indices, coef.indices, drop = FALSE]
1016
1017      if (is.null(estimates))
1018        estimates <- coefvlm(model)[coef.indices]
1019
1020      if (length(refPos) == 1) {
1021        if (length(estimates) != kk)
1022          estimates <- c(0, estimates)
1023        covmat <- rbind(0, cbind(0, covmat))
1024        names(estimates)[1] <-
1025        rownames(covmat)[1] <-
1026        colnames(covmat)[1] <- reference.name
1027        if (refPos != 1) {
1028          perm <- if (refPos == kk) c(2:kk, 1) else
1029                  c(2:refPos, 1, (refPos + 1):kk)
1030          estimates <- estimates[perm]
1031          covmat <- covmat[perm, perm, drop = FALSE]
1032        }
1033      }
1034    }  # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
1035
1036
1037    return(Recall(covmat,
1038                  factorname = factorname,
1039                  which.linpred = which.linpred,
1040                  coef.indices = coef.indices.saved,
1041                  labels = labels,
1042                  dispersion = dispersion,
1043                  estimates = estimates
1044                  )
1045           )
1046  } else { # ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1047
1048    covmat <- object
1049    if (length(labels))
1050      rownames(covmat) <- colnames(covmat) <- labels
1051    if ((LLL <- dim(covmat)[1]) <= 2)
1052      stop("This function works only for factors with 3 ",
1053           "or more levels")
1054  }  # ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1055
1056
1057
1058
1059  allvcov <- covmat
1060  for (ilocal in 1:LLL)
1061    for (jlocal in ilocal:LLL)
1062      allvcov[ilocal, jlocal] <-
1063      allvcov[jlocal, ilocal] <- covmat[ilocal, ilocal] +
1064                                 covmat[jlocal, jlocal] -
1065                                 covmat[ilocal, jlocal] * 2
1066
1067  diag(allvcov) <- rep_len(1.0, LLL)  # Any positive value should do
1068
1069
1070  wmat   <- matrix(1.0, LLL, LLL)
1071  diag(wmat) <- sqrt( .Machine$double.eps )
1072
1073
1074  logAllvcov <- log(allvcov)
1075  attr(logAllvcov, "Prior.Weights") <- wmat
1076  attr(logAllvcov, "estimates")     <- estimates
1077  attr(logAllvcov, "coef.indices")  <- coef.indices
1078  attr(logAllvcov, "factorname")    <- factorname
1079  attr(logAllvcov, "regularVar")    <- diag(covmat)
1080  attr(logAllvcov, "which.linpred") <- which.linpred
1081
1082  logAllvcov
1083}  # End of Qvar()
1084
1085
1086
1087
1088
1089
1090
1091
1092WorstErrors <- function(qv.object) {
1093  stop("20110729; does not work")
1094
1095  reducedForm <- function(covmat, qvmat) {
1096    nlevels <- dim(covmat)[1]
1097    firstRow <- covmat[1, ]
1098    ones <- rep_len(1, nlevels)
1099    J <- outer(ones, ones)
1100    notzero <- 2:nlevels
1101    r.covmat <- covmat + (firstRow[1]*J) -
1102                         outer(firstRow, ones) -
1103                         outer(ones, firstRow)
1104    r.covmat <- r.covmat[notzero, notzero]
1105    qv1 <- qvmat[1, 1]
1106    r.qvmat <- (qvmat + qv1*J)[notzero, notzero]
1107    list(r.covmat, r.qvmat)
1108  }
1109  covmat <- qv.object$covmat
1110  qvmat <- diag(qv.object$qvframe$quasiVar)
1111  r.form <- reducedForm(covmat, qvmat)
1112  r.covmat <- r.form[[1]]
1113  r.qvmat <- r.form[[2]]
1114  inverse.sqrt <- solve(chol(r.covmat))
1115  evalues <- eigen(t(inverse.sqrt) %*% r.qvmat %*% inverse.sqrt,
1116                   symmetric = TRUE)$values
1117  sqrt(c(min(evalues), max(evalues))) - 1
1118}
1119
1120
1121
1122
1123IndentPrint <- function(object, indent = 4, ...) {
1124  stop("20110729; does not work")
1125
1126  zz <- ""
1127  tc <- textConnection("zz", "w", local = TRUE)
1128  sink(tc)
1129  try(print(object, ...))
1130  sink()
1131  close(tc)
1132  indent <- paste(rep_len(" ", indent), sep = "", collapse = "")
1133  cat(paste(indent, zz, sep = ""), sep = "\n")
1134}
1135
1136
1137
1138Print.qv <- function(x, ...) {
1139  stop("20110729; does not work")
1140
1141}
1142
1143
1144
1145
1146summary.qvar <- function(object, ...) {
1147
1148
1149  relerrs <- 1 - sqrt(exp(residuals(object, type = "response")))
1150  diag(relerrs) <- NA
1151
1152  minErrSimple <- round(100 * min(relerrs, na.rm = TRUE), 1)
1153  maxErrSimple <- round(100 * max(relerrs, na.rm = TRUE), 1)
1154
1155
1156
1157  estimates <- c(object@extra$attributes.y$estimates)
1158  if (!length(names(estimates)) &&
1159      is.matrix(object@extra$attributes.y$estimates))
1160    names( estimates) <- rownames(object@extra$attributes.y$estimates)
1161  if (!length(names(estimates)))
1162    names( estimates) <- param.names("Level", length(estimates))
1163
1164
1165  regularVar <- c(object@extra$attributes.y$regularVar)
1166  QuasiVar <- exp(diag(fitted(object))) / 2
1167  QuasiSE  <- sqrt(QuasiVar)
1168
1169
1170  structure(list(estimate = estimates,
1171                 SE            = sqrt(regularVar),
1172                 minErrSimple  = minErrSimple,
1173                 maxErrSimple  = maxErrSimple,
1174                 quasiSE  = QuasiSE,
1175                 object   = object,
1176                 quasiVar = QuasiVar),
1177            class = "summary.qvar")
1178}
1179
1180
1181
1182
1183print.summary.qvar <- function(x, ...) {
1184
1185  object <- x$object
1186  minErrSimple  <- x$minErrSimple
1187  maxErrSimple  <- x$maxErrSimple
1188
1189  x$minErrSimple <- NULL
1190  x$maxErrSimple <- NULL
1191  x$object <- NULL
1192
1193
1194  if (length(cl <- object@call)) {
1195      cat("Call:\n")
1196      dput(cl)
1197  }
1198
1199
1200  facname <- c(object@extra$attributes.y$factorname)
1201  if (length(facname))
1202    cat("Factor name: ", facname, "\n")
1203
1204
1205  if (length(object@dispersion))
1206    cat("\nDispersion: ", object@dispersion, "\n\n")
1207
1208  x <- as.data.frame(c(x))
1209  print.data.frame(x)
1210
1211
1212    cat("\nWorst relative errors in SEs of simple contrasts (%): ",
1213        minErrSimple, ", ", maxErrSimple, "\n")
1214
1215  invisible(x)
1216}
1217
1218
1219
1220
1221qvar <- function(object, se = FALSE, ...) {
1222
1223
1224
1225  if (!inherits(object, "rcim") && !inherits(object, "rcim0"))
1226    warning("argument 'object' should be an 'rcim' or 'rcim0' object. ",
1227            "This call may fail.")
1228
1229  if (!(object@family@vfamily %in% c("uninormal", "normal1")))
1230    warning("argument 'object' does not seem to have used ",
1231            "a 'uninormal' family.")
1232
1233  if (!any(object@misc$link == "explink"))
1234    warning("argument 'object' does not seem to have used ",
1235            "a 'explink' link function.")
1236
1237  quasiVar <- diag(predict(object)[, c(TRUE, FALSE)]) / 2
1238  if (se) sqrt(quasiVar) else quasiVar
1239}
1240
1241
1242
1243
1244
1245
1246
1247
1248plotqvar <-
1249qvplot   <-  function(object,
1250                      interval.width = 2,
1251                      ylab = "Estimate",
1252                      xlab = NULL,  # x$factorname,
1253                      ylim = NULL,
1254                      main = "",
1255                      level.names = NULL,
1256                      conf.level = 0.95,
1257                      warn.ratio = 10,
1258                      border = "transparent",  # None
1259                      points.arg = TRUE,
1260                      length.arrows = 0.25, angle = 30,
1261                      lwd = par()$lwd,
1262                      scol = par()$col,
1263                      slwd = par()$lwd,
1264                      slty = par()$lty,
1265                      ...) {
1266
1267
1268
1269
1270  if (!is.numeric(interval.width) &&
1271      !is.numeric(conf.level))
1272    stop("at least one of arguments 'interval.width' and 'conf.level' ",
1273          "should be numeric")
1274
1275
1276
1277
1278
1279  if (!any("uninormal" %in% object@family@vfamily))
1280    stop("argument 'object' dos not appear to be a ",
1281         "rcim(, uninormal) object")
1282
1283  estimates <- c(object@extra$attributes.y$estimates)
1284  if (!length(names(estimates)) &&
1285      is.matrix(object@extra$attributes.y$estimates))
1286    names(estimates) <- rownames(object@extra$attributes.y$estimates)
1287
1288
1289  if (length(level.names) == length(estimates)) {
1290    names(estimates) <- level.names
1291  } else if (!length(names(estimates)))
1292    names(estimates) <- param.names("Level", length(estimates))
1293
1294
1295
1296
1297  QuasiVar <- exp(diag(fitted(object))) / 2
1298  QuasiSE  <- sqrt(QuasiVar)
1299
1300  if (!is.numeric(estimates))
1301    stop("Cannot plot, because there are no 'proper' ",
1302          "parameter estimates")
1303  if (!is.numeric(QuasiSE))
1304    stop("Cannot plot, because there are no ",
1305         "quasi standard errors")
1306
1307
1308
1309  faclevels <- factor(names(estimates), levels = names(estimates))
1310
1311
1312  xvalues <- seq(along = faclevels)
1313  tops  <- estimates + interval.width * QuasiSE
1314  tails <- estimates - interval.width * QuasiSE
1315
1316
1317
1318
1319  if (is.numeric(conf.level)) {
1320    zedd <- abs(qnorm((1 - conf.level) / 2))
1321    lsd.tops  <- estimates + zedd * QuasiSE / sqrt(2)
1322    lsd.tails <- estimates - zedd * QuasiSE / sqrt(2)
1323    if (max(QuasiSE) / min(QuasiSE) > warn.ratio)
1324      warning("Quasi SEs appear to be quite different... the ",
1325              "LSD intervals may not be very accurate")
1326  } else {
1327    lsd.tops  <- NULL
1328    lsd.tails <- NULL
1329  }
1330
1331
1332
1333
1334  if (is.null(ylim))
1335    ylim <- range(c(tails, tops, lsd.tails, lsd.tops),
1336                  na.rm = TRUE)
1337
1338  if (is.null(xlab))
1339    xlab <- "Factor level"
1340
1341  plot(faclevels, estimates,
1342       border = border,
1343       ylim = ylim, xlab = xlab, ylab = ylab,
1344       lwd = lwd,
1345       main = main, ...)
1346
1347
1348  if (points.arg)
1349    points(estimates, ...)
1350
1351
1352  if (is.numeric(interval.width)) {
1353    segments(xvalues, tails, xvalues, tops,
1354             col = scol, lty = slty, lwd = slwd)
1355  }
1356
1357
1358  if (is.numeric(conf.level)) {
1359    arrows(xvalues, lsd.tails, xvalues, lsd.tops,
1360           col = scol, lty = slty, lwd = slwd, code = 3,
1361           length = length.arrows, angle = angle)
1362
1363  }
1364
1365
1366
1367
1368  if (any(slotNames(object) == "post")) {
1369    object@post$estimates  <- estimates
1370    object@post$xvalues    <- xvalues
1371    if (is.numeric(interval.width)) {
1372      object@post$tails <- tails
1373      object@post$tops  <- tops
1374    }
1375    if (is.numeric(conf.level)) {
1376      object@post$lsd.tails <- lsd.tails
1377      object@post$lsd.tops  <- lsd.tops
1378    }
1379  }
1380
1381
1382  invisible(object)
1383}
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398