1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17profilevglm <-
18  function(object, which = 1:p.vlm, alpha = 0.01,
19           maxsteps = 10, del = zmax/5, trace = NULL, ...) {
20
21
22
23
24  Pnames <- names(B0 <- coef(object))
25  nonA <- !is.na(B0)
26  if (any(is.na(B0)))
27    stop("currently cannot handle NA-valued regression coefficients")
28  pv0 <- t(as.matrix(B0))  # 1 x p.vlm
29
30
31
32  p.vlm <- length(Pnames)
33  if (is.character(which))
34    which <- match(which, Pnames)
35  summ <- summary(object)
36  std.err <- coef(summ)[, "Std. Error", drop = FALSE]
37
38
39
40  M <- npred(object)
41  Xm2 <- model.matrix(object, type = "lm2")  # Could be a 0 x 0 matrix
42  if (!length(Xm2))
43     Xm2 <- NULL  # Make sure. This is safer
44
45
46
47
48  clist <- constraints(object, type = "term")  # type = c("lm", "term")
49
50
51
52
53  mf <- model.frame(object)
54
55  Y <- model.response(mf)
56  if (!is.factor(Y))
57    Y <- as.matrix(Y)
58
59
60  n.lm <- nobs(object, type = "lm")
61  OOO <- object@offset
62  if (!length(OOO) || all(OOO == 0))
63    OOO <- matrix(0, n.lm, M)
64
65
66
67  mt <- attr(mf, "terms")
68
69
70
71
72  Wts <- model.weights(mf)
73  if (length(Wts) == 0L)
74    Wts <- rep(1, n.lm)  # Safest (uses recycling and is a vector)
75  Original.de <- deviance(object)  # Could be NULL
76  if (!(use.de <- is.Numeric(Original.de)))
77    Original.ll <- logLik(object)
78  DispersionParameter <- summ@dispersion
79  if (!all(DispersionParameter == 1))
80    stop("Currently can only handle dispersion parameters ",
81         "that are equal to 1")
82  X.lm  <- model.matrix(object, type =  "lm")
83  X.vlm <- model.matrix(object, type = "vlm")
84  fam <- object@family
85
86
87
88
89
90  quasi.type <- if (length(tmp3 <- fam@infos()$quasi.type))
91    tmp3 else FALSE
92  if (quasi.type)
93    stop("currently this function cannot handle quasi-type models",
94         " or models with an estimated dispersion parameter")
95
96
97  zmax <- sqrt(qchisq(1 - alpha, 1))
98  profName <- "z"
99
100
101  prof <- vector("list", length = length(which))
102  names(prof) <- Pnames[which]
103
104
105  for (i in which) {
106    zi <- 0
107    pvi <- pv0
108    aa <- nonA
109    aa[i] <- FALSE
110    X.vlm.i <- X.vlm[, aa, drop = FALSE]
111    X.lm.i  <-  X.lm  # Try this
112
113
114 # This is needed by vglm.fit():
115    attr(X.vlm.i, "assign") <- attr(X.vlm, "assign")  # zz; this is wrong!
116    attr( X.lm.i, "assign") <- attr( X.lm, "assign")
117
118
119    if (is.logical(trace))
120      object@control$trace <- trace
121
122
123    pnamesi <- Pnames[i]
124    for (sgn in c(-1, 1)) {
125      if (is.logical(trace) && trace)
126        message("\nParameter: ", pnamesi, " ",
127                c("down", "up")[(sgn + 1)/2 + 1])
128      step <- 0
129      zedd <- 0
130      LPmat <- matrix(c(X.vlm[, nonA, drop = FALSE] %*% B0[nonA]),
131                      n.lm, M, byrow = TRUE) + OOO
132
133
134      while ((step <- step + 1) < maxsteps &&
135             abs(zedd) < zmax) {
136        betai <- B0[i] + sgn * step * del * std.err[Pnames[i], 1]
137        ooo <- OOO + matrix(X.vlm[, i] * betai, n.lm, M, byrow = TRUE)
138
139
140
141
142        fm <- vglm.fit(x = X.lm.i,  # Possibly use X.lm.i or else X.lm
143                       y = Y, w = Wts,
144                       X.vlm.arg = X.vlm.i,  # X.vlm,
145                       Xm2 = Xm2, Terms = mt,
146                       constraints = clist, extra = object@extra,
147                       etastart = LPmat,
148                       offset = ooo, family = fam,
149                       control = object@control)
150
151
152
153        fmc <- fm$coefficients
154        LPmat <- matrix(X.vlm.i %*% fmc, n.lm, M, byrow = TRUE) + ooo
155        ri <- pv0
156        ri[, names(fmc)] <- fmc  # coef(fm)
157        ri[, pnamesi] <- betai
158        pvi <- rbind(pvi, ri)
159        zee <- if (use.de) {
160          fm$crit.list[["deviance"]] - Original.de
161        } else {
162          2 * (Original.ll - fm$crit.list[["loglikelihood"]])
163        }
164        if (zee > -1e-3) {
165          zee <- max(zee, 0)
166        } else {
167          stop("profiling has found a better solution, ",
168               "so original fit had not converged")
169        }
170        zedd <- sgn * sqrt(zee)
171        zi <- c(zi, zedd)
172      }  # while
173    }  # for sgn
174    si. <- order(zi)
175    prof[[pnamesi]] <- structure(data.frame(zi[si.]), names = profName)
176    prof[[pnamesi]]$par.vals <- pvi[si., ,drop = FALSE]
177  }  # for i
178
179
180  val <- structure(prof, original.fit = object, summary = summ)
181  class(val) <- c("profile.glm", "profile")
182  val
183}
184
185
186
187
188
189
190if (!isGeneric("profile"))
191    setGeneric("profile",
192               function(fitted, ...)
193               standardGeneric("profile"),
194           package = "VGAM")
195
196
197setMethod("profile", "vglm",
198          function(fitted, ...)
199          profilevglm(object = fitted, ...))
200
201
202
203
204
205
206
207vplot.profile <-
208  function(x, ...) {
209  nulls <- sapply(x, is.null)
210  if (all(nulls)) return(NULL)
211  x <- x[!nulls]
212  nm <- names(x)
213  nr <- ceiling(sqrt(length(nm)))
214  oldpar <- par(mfrow = c(nr, nr))
215  on.exit(par(oldpar))
216  for (nm in names(x)) {
217    tau <- x[[nm]][[1L]]
218    parval <- x[[nm]][[2L]][, nm]
219    dev.hold()
220    plot(parval, tau, xlab = nm, ylab = "tau", type = "n")
221    if (sum(tau == 0) == 1) points(parval[tau == 0], 0, pch = 3)
222    splineVals <- spline(parval, tau)
223    lines(splineVals$x, splineVals$y)
224    dev.flush()
225  }
226}
227
228
229
230vpairs.profile <-
231function(x, colours = 2:3, ...) {
232  parvals <- lapply(x, "[[", "par.vals")
233
234
235  rng <- apply(do.call("rbind", parvals), 2L, range, na.rm = TRUE)
236  Pnames <- colnames(rng)
237  npar <- length(Pnames)
238  coefs <- coef(attr(x, "original.fit"))
239  form <- paste(as.character(formula(attr(x, "original.fit")))[c(2, 1, 3)],
240                collapse = "")
241  oldpar <- par(mar = c(0, 0, 0, 0), mfrow = c(1, 1),
242                oma = c(3, 3, 6, 3), las = 1)
243  on.exit(par(oldpar))
244  fin <- par("fin")
245  dif <- (fin[2L] - fin[1L])/2
246  adj <- if (dif > 0) c(dif, 0, dif, 0) else c(0, -dif, 0, -dif)
247  par(omi = par("omi") + adj)
248  cex <- 1 + 1/npar
249  frame()
250  mtext(form, side = 3, line = 3, cex = 1.5, outer = TRUE)
251  del <- 1/npar
252  for (i in 1L:npar) {
253    ci <- npar - i
254    pi <- Pnames[i]
255    for (j in 1L:npar) {
256      dev.hold()
257      pj <- Pnames[j]
258      par(fig = del * c(j - 1, j, ci, ci + 1))
259      if (i == j) {
260        par(new = TRUE)
261        plot(rng[, pj], rng[, pi], axes = FALSE,
262             xlab = "", ylab = "", type = "n")
263        op <- par(usr = c(-1, 1, -1, 1))
264        text(0, 0, pi, cex = cex, adj = 0.5)
265        par(op)
266      } else {
267        col <- colours
268        if (i < j) col <- col[2:1]
269        if (!is.null(parvals[[pj]])) {
270          par(new = TRUE)
271          plot(spline(x <- parvals[[pj]][, pj],
272                      y <- parvals[[pj]][, pi]),
273               type = "l", xlim = rng[, pj],
274               ylim = rng[, pi], axes = FALSE,
275               xlab = "", ylab = "", col = col[2L])
276          pu <- par("usr")
277          smidge <- 2/100 * (pu[4L] - pu[3L])
278          segments(x, pmax(pu[3L], y - smidge),
279                   x, pmin(pu[4L], y + smidge))
280        } else
281          plot(rng[, pj], rng[, pi], axes = FALSE,
282               xlab = "", ylab = "", type = "n")
283        if (!is.null(parvals[[pi]])) {
284          lines(x <- parvals[[pi]][, pj],
285                y <- parvals[[pi]][, pi],
286                type = "l", col = col[1L])
287          pu <- par("usr")
288          smidge <- 2/100 * (pu[2L] - pu[1L])
289          segments(pmax(pu[1L], x - smidge), y,
290                   pmin(pu[2L], x + smidge), y)
291        }
292        points(coefs[pj], coefs[pi], pch = 3, cex = 3)
293      }
294      if (i == npar) axis(1)
295      if (j == 1) axis(2)
296      if (i == 1) axis(3)
297      if (j == npar) axis(4)
298      dev.flush()
299    }
300  }
301  par(fig = c(0, 1, 0, 1))
302  invisible(x)
303}
304
305
306
307
308
309