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 18N.hat.posbernoulli <- 19 function(eta, link, earg = list(), 20 R = NULL, w = NULL, 21 X.vlm = NULL, Hlist = NULL, 22 extra = list(), 23 model.type = c("0", "b", "t", "tb") 24 ) { 25 26 27 28 29 30 if (!is.null(w) && !all(w[1] == w)) 31 warning("estimate of N may be wrong when prior weights ", 32 "are not all the same") 33 34 model.type <- match.arg(model.type, c("0", "b", "t", "tb"))[1] 35 if (!is.matrix(eta)) 36 eta <- as.matrix(eta) # May be needed for "0" 37 38 tau <- 39 switch(model.type, 40 "0" = extra$tau, 41 "b" = extra$tau, 42 "t" = ncol(eta), 43 "tb" = (ncol(eta) + 1) / 2) 44 if (length(extra$tau) && extra$tau != tau) 45 warning("variable 'tau' is mistaken") # Checking only 46 47 48 jay.index <- 49 switch(model.type, 50 "0" = rep_len(1, tau), 51 "b" = rep_len(1, tau), # Subset: 2 out of 1:2 52 "t" = 1:tau, # All of them 53 "tb" = 1:tau) # Subset: first tau of them out of M = 2*tau-2 54 55 prc <- eta2theta(eta[, jay.index], link, earg = earg) # cap.probs 56 prc <- as.matrix(prc) # Might be needed for Mtb(tau=2). 57 58 59 60 if (FALSE && model.type == "tb") { 61 if (tau == 2) 62 prc <- cbind(prc, 1 - prc) 63 if (tau > 3) 64 stop("cannot handle tau > 3 yet") 65 jay.index <- 1:tau # 'Restore' it coz its used below. zz?? 66 } 67 68 QQQ <- exp(rowSums(log1p(-prc))) 69 pibbeta <- exp(log1p(-QQQ)) # One.minus.QQQ 70 N.hat <- sum(1 / pibbeta) # Point estimate 71 ss2 <- sum(QQQ / pibbeta^2) # Assumes bbeta is known 72 73 74 if (length(extra$p.small) && 75 any(pibbeta < extra$p.small) && 76 !extra$no.warning) 77 warning("The abundance estimation for this model can be unstable") 78 79 80 if (length(R)) { 81 82 dvect <- matrix(0, length(pibbeta), ncol = ncol(X.vlm)) 83 M <- nrow(Hlist[[1]]) 84 n.lm <- nrow(X.vlm) / M # Number of rows of the LM matrix 85 dprc.deta <- dtheta.deta(prc, link, earg = earg) 86 Hmatrices <- matrix(c(unlist(Hlist)), nrow = M) 87 for (jay in 1:tau) { 88 linpred.index <- jay.index[jay] 89 Index0 <- Hmatrices[linpred.index, ] != 0 90 X.lm.jay <- X.vlm[(0:(n.lm - 1)) * M + linpred.index, 91 Index0, 92 drop = FALSE] 93 94 dvect[, Index0] <- 95 dvect[, Index0] + 96 (QQQ / (1-prc[, jay])) * dprc.deta[, jay] * X.lm.jay 97 } 98 99 100 dvect <- dvect * (-1 / pibbeta^2) 101 dvect <- colSums(dvect) # Now a vector 102 103 ncol.X.vlm <- nrow(R) 104 rinv <- diag(ncol.X.vlm) 105 rinv <- backsolve(R, rinv) 106 rowlen <- drop(((rinv^2) %*% rep_len(1, ncol.X.vlm))^0.5) 107 covun <- rinv %*% t(rinv) 108 vecTF <- FALSE 109 for (jay in 1:tau) { 110 linpred.index <- jay.index[jay] 111 vecTF <- vecTF | (Hmatrices[linpred.index, ] != 0) 112 } 113 vecTF.index <- (seq_along(vecTF))[vecTF] 114 covun <- covun[vecTF.index, vecTF.index, drop = FALSE] 115 dvect <- dvect[vecTF.index, drop = FALSE] 116 } 117 118 list(N.hat = N.hat, 119 SE.N.hat = if (length(R)) 120 c(sqrt(ss2 + t(dvect) %*% covun %*% dvect)) else 121 c(sqrt(ss2)) 122 ) 123} # N.hat.posbernoulli 124 125 126 127 128 129 130 aux.posbernoulli.t <- 131 function(y, check.y = FALSE, 132 rename = TRUE, 133 name = "bei") { 134 135 136 137 138 139 140 y <- as.matrix(y) 141 if ((tau <- ncol(y)) == 1) 142 stop("argument 'y' needs to be a matrix with at least two columns") 143 if (check.y) { 144 if (!all(y == 0 | y == 1 | y == 1/tau | is.na(y))) 145 stop("response 'y' must contain 0s and 1s only") 146 } 147 148 zeddij <- cbind(0, t(apply(y, 1, cumsum))) # tau + 1 columns 149 zij <- (0 + (zeddij > 0))[, 1:tau] # 0 or 1. 150 if (rename) { 151 colnames(zij) <- param.names(name, ncol(y)) 152 } else { 153 if (length(colnames(y))) 154 colnames(zij) <- colnames(y) 155 } 156 157 158 cp1 <- numeric(nrow(y)) 159 for (jay in tau:1) 160 cp1[y[, jay] > 0] <- jay 161 if (any(cp1 == 0)) 162 warning("some individuals were never captured!") 163 164 yr1i <- zeddij[, tau + 1] - 1 165 list(cap.hist1 = zij, # A matrix of the same dimension as 'y' 166 cap1 = cp1, # Aka ti1 167 y0i = cp1 - 1, 168 yr0i = tau - cp1 - yr1i, 169 yr1i = yr1i) 170} # aux.posbernoulli.t 171 172 173 174 175 176 177 178 179 180rposbern <- 181 function(n, nTimePts = 5, pvars = length(xcoeff), 182 xcoeff = c(-2, 1, 2), 183 Xmatrix = NULL, # If is.null(Xmatrix) then it is created 184 cap.effect = 1, 185 is.popn = FALSE, 186 link = "logitlink", 187 earg.link = FALSE) { 188 189 190 191 192 193 194 195 196 197 use.n <- if ((length.n <- length(n)) > 1) length.n else 198 if (!is.Numeric(n, integer.valued = TRUE, 199 length.arg = 1, positive = TRUE)) 200 stop("bad input for argument 'n'") else n 201 orig.n <- use.n 202 if (!is.popn) 203 use.n <- 1.50 * use.n + 100 # Bigger due to rejections 204 205 if (pvars == 0) 206 stop("argument 'pvars' must be at least one") 207 if (pvars > length(xcoeff)) 208 stop("argument 'pvars' is too high") 209 210 211 if (earg.link) { 212 earg <- link 213 } else { 214 link <- as.list(substitute(link)) 215 earg <- link2list(link) 216 } 217 link <- attr(earg, "function.name") 218 219 220 cap.effect.orig <- cap.effect 221 222 223 Ymatrix <- matrix(0, use.n, nTimePts, 224 dimnames = list(as.character(1:use.n), 225 param.names("y", nTimePts))) 226 227 CHmatrix <- matrix(0, use.n, nTimePts, 228 dimnames = list(as.character(1:use.n), 229 param.names("ch", nTimePts))) 230 231 232 if (is.null(Xmatrix)) { 233 Xmatrix <- cbind(x1 = rep_len(1.0, use.n)) 234 if (pvars > 1) 235 Xmatrix <- cbind(Xmatrix, 236 matrix(runif(n = use.n * (pvars-1)), 237 use.n, pvars - 1, 238 dimnames = list(as.character(1:use.n), 239 param.names("x", pvars)[-1]))) 240 } 241 242 243 lin.pred.baseline <- xcoeff[1] 244 if (pvars > 1) 245 lin.pred.baseline <- lin.pred.baseline + 246 Xmatrix[, 2:pvars, drop = FALSE] %*% 247 xcoeff[2:pvars] 248 sumrowy <- rep_len(0, use.n) 249 cap.effect <- rep_len(cap.effect.orig, use.n) 250 251 for (jlocal in 1:nTimePts) { 252 CHmatrix[, jlocal] <- as.numeric(sumrowy > 0) 253 254 caught.before.TF <- (CHmatrix[, jlocal] > 0) 255 lin.pred <- lin.pred.baseline + caught.before.TF * cap.effect 256 257 Ymatrix[, jlocal] <- 258 rbinom(use.n, size = 1, 259 prob = eta2theta(lin.pred, link = link, earg = earg)) 260 261 sumrowy <- sumrowy + Ymatrix[, jlocal] 262 } 263 264 265 266 index0 <- (sumrowy == 0) 267 if (all(!index0)) 268 stop("bug in this code: cannot handle no animals being caught") 269 Ymatrix <- Ymatrix[!index0, , drop = FALSE] 270 Xmatrix <- Xmatrix[!index0, , drop = FALSE] 271 CHmatrix <- CHmatrix[!index0, , drop = FALSE] 272 273 274 275 276 ans <- data.frame(Ymatrix, Xmatrix, CHmatrix # zCHmatrix, 277 ) 278 279 280 if (!is.popn) { 281 ans <- if (nrow(ans) >= orig.n) { 282 ans[1:orig.n, ] 283 } else { 284 rbind(ans, 285 Recall(n = orig.n - nrow(ans), 286 nTimePts = nTimePts, 287 pvars = pvars, 288 xcoeff = xcoeff, 289 cap.effect = cap.effect.orig, 290 link = earg, 291 earg.link = TRUE)) 292 } 293 } 294 295 rownames(ans) <- as.character(1:nrow(ans)) 296 297 attr(ans, "pvars") <- pvars 298 attr(ans, "nTimePts") <- nTimePts 299 attr(ans, "cap.effect") <- cap.effect.orig 300 attr(ans, "is.popn") <- is.popn 301 attr(ans, "n") <- n 302 303 ans 304} # rposbern 305 306 307 308 309 310 311 312dposbern <- function(x, prob, prob0 = prob, log = FALSE) { 313 314 315 x <- as.matrix(x) 316 prob <- as.matrix(prob) 317 prob0 <- as.matrix(prob0) 318 319 if (!is.logical(log.arg <- log) || 320 length(log) != 1) 321 stop("bad input for argument 'log'") 322 rm(log) 323 if (ncol(x) < 2) 324 stop("columns of argument 'x' should be 2 or more") 325 326 327 logAA0 <- rowSums(log1p(-prob0)) 328 AA0 <- exp(logAA0) 329 ell1 <- x * log(prob) + (1 - x) * log1p(-prob) - log1p(-AA0) / ncol(x) 330 if (log.arg) ell1 else exp(ell1) 331} # dposbern 332 333 334 335 336 EIM.posNB.specialp <- 337 function(munb, size, 338 y.max = NULL, # Must be an integer 339 cutoff.prob = 0.995, 340 prob0, df0.dkmat, df02.dkmat2, 341 intercept.only = FALSE, 342 second.deriv = TRUE) { 343 344 345 if (intercept.only) { 346 munb <- munb[1] 347 size <- size[1] 348 prob0 <- prob0[1] 349 df0.dkmat <- df0.dkmat[1] 350 df02.dkmat2 <- df02.dkmat2[1] 351 } 352 353 y.min <- 0 # Same as negbinomial(). A fixed const really 354 355 if (!is.numeric(y.max)) { 356 eff.p <- sort(c(cutoff.prob, 1 - cutoff.prob)) 357 y.max <- max(qgaitnbinom(p = eff.p[2], truncate = 0, 358 size, munb.p = munb)) + 10 359 } 360 361 Y.mat <- if (intercept.only) rbind(y.min:y.max) else 362 matrix(y.min:y.max, length(munb), y.max-y.min+1, 363 byrow = TRUE) 364 neff.row <- ifelse(intercept.only, 1, nrow(Y.mat)) 365 neff.col <- ifelse(intercept.only, length(Y.mat), ncol(Y.mat)) 366 367 if (TRUE) { 368 Y.mat2 <- Y.mat + 1 369 trigg.term0 <- if (intercept.only) { 370 371 372 sum(c(dgaitnbinom(Y.mat2, size[1], munb.p = munb[1], 373 truncate = 0)) * 374 c(trigamma(Y.mat2 + size[1]))) 375 } else { 376 } 377 } # FALSE 378 379 380 381 382 383 384 385 386 trigg.term <- 387 if (TRUE) { 388 answerC <- .C("eimpnbinomspecialp", 389 as.integer(intercept.only), 390 as.double(neff.row), as.double(neff.col), 391 as.double(size), 392 as.double(1 - pgaitnbinom(Y.mat, size, munb.p = munb, truncate = 0 393 )), 394 rowsums = double(neff.row)) 395 answerC$rowsums 396 } 397 398 399 400 mymu <- munb / (1 - prob0) # E(Y) 401 ned2l.dk2 <- trigg.term - munb / (size * (size + munb)) - 402 (mymu - munb) / (munb + size)^2 403 404 if (second.deriv) 405 ned2l.dk2 <- ned2l.dk2 - df02.dkmat2 / (1 - prob0) - 406 (df0.dkmat / (1 - prob0))^2 407 ned2l.dk2 408} # end of EIM.posNB.specialp() 409 410 411 412 413 414 415 416 EIM.posNB.speciald <- 417 function(munb, size, 418 y.min = 1, # 20160201; must be an integer 419 y.max = NULL, # Must be an integer 420 cutoff.prob = 0.995, 421 prob0, df0.dkmat, df02.dkmat2, 422 intercept.only = FALSE, 423 second.deriv = TRUE) { 424 425 426 if (intercept.only) { 427 munb <- munb[1] 428 size <- size[1] 429 prob0 <- prob0[1] 430 df0.dkmat <- df0.dkmat[1] 431 df02.dkmat2 <- df02.dkmat2[1] 432 } 433 434 if (!is.numeric(y.max)) { 435 eff.p <- sort(c(cutoff.prob, 1 - cutoff.prob)) 436 y.max <- max(qgaitnbinom(p = eff.p[2], truncate = 0, 437 size, munb.p = munb)) + 10 438 } 439 440 Y.mat <- if (intercept.only) rbind(y.min:y.max) else 441 matrix(y.min:y.max, length(munb), 442 y.max-y.min+1, byrow = TRUE) 443 trigg.term <- if (intercept.only) { 444 445 sum(c(dgaitnbinom(Y.mat, size[1], munb.p = munb[1], 446 truncate = 0)) * 447 c(trigamma(Y.mat + size[1]))) 448 } else { 449 rowSums(dgaitnbinom(Y.mat, size, munb.p = munb, 450 truncate = 0) * 451 trigamma(Y.mat + size)) 452 } 453 454 mymu <- munb / (1 - prob0) # E(Y) 455 ned2l.dk2 <- trigamma(size) - munb / (size * (size + munb)) - 456 (mymu - munb) / (munb + size)^2 - trigg.term 457 if (second.deriv) 458 ned2l.dk2 <- ned2l.dk2 - df02.dkmat2 / (1 - prob0) - 459 (df0.dkmat / (1 - prob0))^2 460 ned2l.dk2 461} # end of EIM.posNB.speciald() 462 463 464 465 466 467posNBD.Loglikfun2 <- function(munbval, sizeval, 468 y, x, w, extraargs) { 469 sum(c(w) * dgaitnbinom(y, sizeval, munb.p = munbval, 470 truncate = 0, log = TRUE)) 471} 472 473 474 475posnegbinomial.control <- function(save.weights = TRUE, ...) { 476 list(save.weights = save.weights) 477} 478 479 480 481 posnegbinomial <- 482 function( 483 zero = "size", 484 type.fitted = c("mean", "munb", "prob0"), 485 mds.min = 1e-3, 486 nsimEIM = 500, 487 cutoff.prob = 0.999, # higher is better for large 'size' 488 eps.trig = 1e-7, 489 max.support = 4000, # 20160201; I have changed this 490 max.chunk.MB = 30, # max.memory = Inf is allowed 491 lmunb = "loglink", lsize = "loglink", 492 imethod = 1, 493 imunb = NULL, 494 iprobs.y = NULL, # 0.35, 495 gprobs.y = ppoints(8), 496 isize = NULL, 497 gsize.mux = exp(c(-30, -20, -15, -10, -6:3))) { 498 499 500 501 if (!is.Numeric(imethod, length.arg = 1, 502 integer.valued = TRUE, positive = TRUE) || 503 imethod > 2) 504 stop("argument 'imethod' must be 1 or 2") 505 506 507 if (length(isize) && !is.Numeric(isize, positive = TRUE)) 508 stop("bad input for argument 'isize'") 509 510 lmunb <- as.list(substitute(lmunb)) 511 emunb <- link2list(lmunb) 512 lmunb <- attr(emunb, "function.name") 513 514 lsize <- as.list(substitute(lsize)) 515 esize <- link2list(lsize) 516 lsize <- attr(esize, "function.name") 517 518 type.fitted <- match.arg(type.fitted, 519 c("mean", "munb", "prob0"))[1] 520 521 522 if (!is.Numeric(eps.trig, length.arg = 1, 523 positive = TRUE) || eps.trig > 0.001) 524 stop("argument 'eps.trig' must be positive and smaller in value") 525 526 if (!is.Numeric(nsimEIM, length.arg = 1, 527 positive = TRUE, integer.valued = TRUE)) 528 stop("argument 'nsimEIM' must be a positive integer") 529 if (nsimEIM <= 30) 530 warning("argument 'nsimEIM' should be greater than 30, say") 531 532 533 new("vglmff", 534 blurb = c("Positive-negative binomial distribution\n\n", 535 "Links: ", 536 namesof("munb", lmunb, earg = emunb), ", ", 537 namesof("size", lsize, earg = esize), "\n", 538 "Mean: munb / (1 - (size / (size + munb))^size)"), 539 constraints = eval(substitute(expression({ 540 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 541 predictors.names = predictors.names, 542 M1 = 2) 543 }), list( .zero = zero ))), 544 infos = eval(substitute(function(...) { 545 list(M1 = 2, 546 Q1 = 1, 547 expected = TRUE, 548 mds.min = .mds.min , 549 multipleResponses = TRUE, 550 parameters.names = c("munb", "size"), 551 nsimEIM = .nsimEIM , 552 eps.trig = .eps.trig , 553 lmunb = .lmunb , 554 emunb = .emunb , 555 type.fitted = .type.fitted , 556 zero = .zero , 557 lsize = .lsize , 558 esize = .esize ) 559 }, list( .lmunb = lmunb, .lsize = lsize, .isize = isize, 560 .emunb = emunb, .esize = esize, 561 .zero = zero, .nsimEIM = nsimEIM, 562 .eps.trig = eps.trig, 563 .imethod = imethod, 564 .type.fitted = type.fitted, 565 .mds.min = mds.min))), 566 567 initialize = eval(substitute(expression({ 568 M1 <- 2 569 570 temp12 <- 571 w.y.check(w = w, y = y, 572 Is.integer.y = TRUE, 573 Is.positive.y = TRUE, 574 ncol.w.max = Inf, 575 ncol.y.max = Inf, 576 out.wy = TRUE, 577 colsyperw = 1, 578 maximize = TRUE) 579 w <- temp12$w 580 y <- temp12$y 581 582 583 584 M <- M1 * ncol(y) 585 extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species 586 extra$type.fitted <- .type.fitted 587 extra$colnames.y <- colnames(y) 588 589 predictors.names <- c(namesof(param.names("munb", NOS, skip1 = TRUE), 590 .lmunb , earg = .emunb , tag = FALSE), 591 namesof(param.names("size", NOS, skip1 = TRUE), 592 .lsize , earg = .esize , tag = FALSE)) 593 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 594 595 gprobs.y <- .gprobs.y 596 imunb <- .imunb # Default in NULL 597 if (length(imunb)) 598 imunb <- matrix(imunb, n, NOS, byrow = TRUE) 599 600 if (!length(etastart)) { 601 602 603 munb.init <- 604 size.init <- matrix(NA_real_, n, NOS) 605 gprobs.y <- .gprobs.y 606 if (length( .iprobs.y )) 607 gprobs.y <- .iprobs.y 608 gsize.mux <- .gsize.mux # gsize.mux is on a relative scale 609 610 for (jay in 1:NOS) { # For each response 'y_jay'... do: 611 wm.yj <- weighted.mean(y[, jay], w = w[, jay]) 612 munb.init.jay <- if ( .imethod == 1 ) { 613 614 negbinomial.initialize.yj(y[, jay] - 1, 615 w[, jay], gprobs.y = gprobs.y, 616 wm.yj = wm.yj) + 1 - 1/4 617 } else { 618 wm.yj - 1/2 619 } 620 if (length(imunb)) { 621 munb.init.jay <- sample(x = imunb[, jay], 622 size = 10, replace = TRUE) 623 munb.init.jay <- unique(sort(munb.init.jay)) 624 } 625 626 gsize <- gsize.mux * 0.5 * (mean(munb.init.jay) + wm.yj) 627 if (length( .isize )) 628 gsize <- .isize # isize is on an absolute scale 629 630 631 try.this <- 632 grid.search2(munb.init.jay, gsize, 633 objfun = posNBD.Loglikfun2, 634 y = y[, jay], w = w[, jay], 635 ret.objfun = TRUE) # Last value is the loglik 636 munb.init[, jay] <- try.this["Value1"] 637 size.init[, jay] <- try.this["Value2"] 638 } # for (jay ...) 639 640 641 642 643 644 645 646 etastart <- 647 cbind( 648 theta2eta(munb.init , .lmunb , earg = .emunb ), 649 theta2eta(size.init, .lsize , earg = .esize )) 650 etastart <- etastart[, interleave.VGAM(M, M1 = M1), drop = FALSE] 651 } 652 }), list( .lmunb = lmunb, .lsize = lsize, 653 .imunb = imunb, .isize = isize, 654 .emunb = emunb, .esize = esize, 655 .gprobs.y = gprobs.y, .gsize.mux = gsize.mux, 656 .iprobs.y = iprobs.y, 657 .imethod = imethod, 658 .type.fitted = type.fitted ))), 659 660 linkinv = eval(substitute(function(eta, extra = NULL) { 661 NOS <- ncol(eta) / c(M1 = 2) 662 type.fitted <- if (length(extra$type.fitted)) extra$type.fitted else { 663 warning("cannot find 'type.fitted'. ", 664 "Returning the 'mean'.") 665 "mean" 666 } 667 668 type.fitted <- match.arg(type.fitted, 669 c("mean", "munb", "prob0"))[1] 670 671 TF <- c(TRUE, FALSE) 672 munb <- eta2theta(eta[, TF, drop = FALSE], .lmunb , earg = .emunb ) 673 kmat <- eta2theta(eta[, !TF, drop = FALSE], .lsize , earg = .esize ) 674 small.size <- 1e-10 675 if (any(ind4 <- (kmat < small.size))) { 676 warning("estimates of 'size' are very small. ", 677 "Taking evasive action.") 678 kmat[ind4] <- small.size 679 } 680 681 tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb) 682 prob0 <- tempk^kmat 683 oneminusf0 <- 1 - prob0 684 685 smallval <- .mds.min # Something like this is needed 686 687 if (any(big.size <- (munb / kmat < smallval))) { 688 prob0[big.size] <- exp(-munb[big.size]) # The limit as kmat-->Inf 689 oneminusf0[big.size] <- -expm1(-munb[big.size]) 690 } 691 692 ans <- switch(type.fitted, 693 "mean" = munb / oneminusf0, 694 "munb" = munb, 695 "prob0" = prob0) # P(Y=0) 696 label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS) 697 }, list( .lsize = lsize, .lmunb = lmunb, 698 .esize = esize, .emunb = emunb, 699 .mds.min = mds.min ))), 700 last = eval(substitute(expression({ 701 temp0303 <- c(rep_len( .lmunb , NOS), 702 rep_len( .lsize , NOS)) 703 names(temp0303) <- c(param.names("munb", NOS, skip1 = TRUE), 704 param.names("size", NOS, skip1 = TRUE)) 705 temp0303 <- temp0303[interleave.VGAM(M, M1 = M1)] 706 misc$link <- temp0303 # Already named 707 708 misc$earg <- vector("list", M1*NOS) 709 names(misc$earg) <- names(misc$link) 710 for (ii in 1:NOS) { 711 misc$earg[[M1*ii-1]] <- .emunb 712 misc$earg[[M1*ii ]] <- .esize 713 } 714 715 misc$max.chunk.MB <- .max.chunk.MB 716 misc$cutoff.prob <- .cutoff.prob 717 misc$imethod <- .imethod 718 misc$nsimEIM <- .nsimEIM 719 misc$expected <- TRUE 720 misc$multipleResponses <- TRUE 721 }), list( .lmunb = lmunb, .lsize = lsize, 722 .emunb = emunb, .esize = esize, 723 .cutoff.prob = cutoff.prob, 724 .max.chunk.MB = max.chunk.MB, 725 .nsimEIM = nsimEIM, .imethod = imethod ))), 726 loglikelihood = eval(substitute( 727 function(mu, y, w, residuals = FALSE, eta, 728 extra = NULL, 729 summation = TRUE) { 730 TFvec <- c(TRUE, FALSE) 731 munb <- eta2theta(eta[, TFvec, drop = FALSE], .lmunb , earg = .emunb ) 732 kmat <- eta2theta(eta[, !TFvec, drop = FALSE], .lsize , earg = .esize ) 733 if (residuals) { 734 stop("loglikelihood residuals not implemented yet") 735 } else { 736 ll.elts <- 737 c(w) * dgaitnbinom(y, kmat, munb.p = munb, 738 truncate = 0, log = TRUE) 739 if (summation) { 740 sum(ll.elts) 741 } else { 742 ll.elts 743 } 744 } 745 }, list( .lmunb = lmunb, .lsize = lsize, 746 .emunb = emunb, .esize = esize ))), 747 748 vfamily = c("posnegbinomial"), 749 750 751 752 753 simslot = eval(substitute( 754 function(object, nsim) { 755 756 pwts <- if (length(pwts <- object@prior.weights) > 0) 757 pwts else weights(object, type = "prior") 758 if (any(pwts != 1)) 759 warning("ignoring prior weights") 760 eta <- predict(object) 761 munb <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], 762 .lmunb, earg = .emunb ) 763 kmat <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], 764 .lsize, earg = .esize ) 765 rgaitnbinom(nsim * length(munb), kmat, munb.p = munb, truncate = 0) 766 }, list( .lmunb = lmunb, .lsize = lsize, 767 .emunb = emunb, .esize = esize ))), 768 769 770 771 validparams = eval(substitute(function(eta, y, extra = NULL) { 772 munb <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], 773 .lmunb , earg = .emunb ) 774 size <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], 775 .lsize , earg = .esize ) 776 777 small.size.absolute <- 1e-14 # 20160909 778 779 smallval <- .mds.min # .munb.div.size 780 okay1 <- all(is.finite(munb)) && all(0 < munb) && 781 all(is.finite(size)) && all(0 < size) && 782 all(small.size.absolute < size) 783 overdispersion <- if (okay1) all(smallval < munb / size) else FALSE 784 if (!overdispersion) 785 warning("parameter 'size' has very large values relative ", 786 "to 'munb'; ", 787 "try fitting a positive-Poisson ", 788 "model instead.") 789 okay1 && overdispersion 790 }, list( .lmunb = lmunb, .emunb = emunb, 791 .lsize = lsize, .esize = esize, 792 .mds.min = mds.min))), 793 794 795 deriv = eval(substitute(expression({ 796 M1 <- 2 797 NOS <- extra$NOS 798 799 TFvec <- c(TRUE, FALSE) 800 munb <- eta2theta(eta[, TFvec, drop = FALSE], .lmunb , earg = .emunb ) 801 kmat <- eta2theta(eta[, !TFvec, drop = FALSE], .lsize , earg = .esize ) 802 803 804 smallval <- .mds.min # Something like this is needed 805 if (any(big.size <- munb / kmat < smallval)) { 806 if (FALSE) 807 warning("parameter 'size' has very large values; ", 808 "try fitting a positive-Poisson ", 809 "model instead") 810 kmat[big.size] <- munb[big.size] / smallval 811 } 812 813 814 dmunb.deta <- dtheta.deta(munb, .lmunb , earg = .emunb ) 815 dsize.deta <- dtheta.deta(kmat, .lsize , earg = .esize ) 816 817 818 tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb) 819 tempm <- munb / (kmat + munb) 820 prob0 <- tempk^kmat 821 oneminusf0 <- 1 - prob0 822 AA16 <- tempm + log(tempk) 823 df0.dmunb <- -tempk * prob0 824 df0.dkmat <- prob0 * AA16 825 df02.dmunb2 <- prob0 * tempk * (1 + 1/kmat) / (1 + munb/kmat) 826 df02.dkmat2 <- prob0 * ((tempm^2) / kmat + AA16^2) 827 df02.dkmat.dmunb <- -prob0 * (tempm/kmat + AA16) / (1 + munb/kmat) 828 829 830 831 if (any(big.size)) { 832 prob0[big.size] <- exp(-munb[big.size]) # The limit as kmat-->Inf 833 oneminusf0[big.size] <- -expm1(-munb[big.size]) 834 df0.dmunb[big.size] <- -tempk[big.size] * prob0[big.size] 835 df0.dkmat[big.size] <- prob0[big.size] * AA16[big.size] 836 df02.dmunb2[big.size] <- prob0[big.size] * tempk[big.size] * 837 (1 + 1/kmat[big.size]) / (1 + smallval) 838 df02.dkmat2[big.size] <- prob0[big.size] * 839 ((tempm[big.size])^2 / kmat[big.size] + AA16[big.size]^2) 840 df02.dkmat.dmunb[big.size] <- -prob0[big.size] * 841 (tempm[big.size]/kmat[big.size] + AA16[big.size]) / (1 + smallval) 842 } 843 844 845 846 847 smallno <- 1e-6 848 if (TRUE && any(near.boundary <- oneminusf0 < smallno)) { 849 warning("solution near the boundary; either there is no need ", 850 "to fit a positive NBD or the distribution is centred ", 851 "on the value 1") 852 oneminusf0[near.boundary] <- smallno 853 prob0[near.boundary] <- 1 - oneminusf0[near.boundary] 854 } 855 856 857 858 859 dl.dmunb <- y / munb - (1 + y/kmat) / (1 + munb/kmat) + 860 df0.dmunb / oneminusf0 861 dl.dsize <- digamma(y + kmat) - digamma(kmat) - 862 (y - munb) / (munb + kmat) + log(tempk) + 863 df0.dkmat / oneminusf0 864 865 866 if (any(big.size)) { 867 } 868 869 870 871 myderiv <- c(w) * cbind(dl.dmunb * dmunb.deta, 872 dl.dsize * dsize.deta) 873 myderiv[, interleave.VGAM(M, M1 = M1)] 874 }), list( .lmunb = lmunb, .lsize = lsize, 875 .emunb = emunb, .esize = esize, 876 .mds.min = mds.min ))), 877 878 879 weight = eval(substitute(expression({ 880 wz <- matrix(0, n, M+M-1) 881 mymu <- munb / oneminusf0 # Is the same as 'mu', == E(Y) 882 883 max.support <- .max.support 884 max.chunk.MB <- .max.chunk.MB 885 886 887 888 889 890 ind2 <- matrix(FALSE, n, NOS) # Used for SFS 891 for (jay in 1:NOS) { 892 eff.p <- sort(c( .cutoff.prob , 1 - .cutoff.prob )) 893 Q.mins <- 1 894 Q.maxs <- qgaitnbinom(p = eff.p[2], truncate = 0, 895 kmat[, jay], munb.p = munb[, jay]) + 10 896 897 898 eps.trig <- .eps.trig 899 Q.MAXS <- pmax(10, ceiling(1 / sqrt(eps.trig))) # 900 Q.maxs <- pmin(Q.maxs, Q.MAXS) 901 902 903 ind1 <- if (max.chunk.MB > 0) 904 (Q.maxs - Q.mins < max.support) else FALSE 905 if ((NN <- sum(ind1)) > 0) { 906 Object.Size <- NN * 8 * max(Q.maxs - Q.mins) / (2^20) 907 n.chunks <- if (intercept.only) 1 else 908 max(1, ceiling( Object.Size / max.chunk.MB)) 909 chunk.rows <- ceiling(NN / n.chunks) 910 ind2[, jay] <- ind1 # Save this 911 wind2 <- which(ind1) 912 913 914 upr.ptr <- 0 915 lwr.ptr <- upr.ptr + 1 916 while (lwr.ptr <= NN) { 917 upr.ptr <- min(upr.ptr + chunk.rows, NN) 918 sind2 <- wind2[lwr.ptr:upr.ptr] 919 wz[sind2, M1*jay] <- 920 EIM.posNB.specialp(munb = munb[sind2, jay], 921 size = kmat[sind2, jay], 922 y.max = max(Q.maxs[sind2]), 923 cutoff.prob = .cutoff.prob , 924 prob0 = prob0[sind2, jay], 925 df0.dkmat = df0.dkmat[sind2, jay], 926 df02.dkmat2 = df02.dkmat2[sind2, jay], 927 intercept.only = intercept.only) 928 929 930 if (any(eim.kk.TF <- wz[sind2, M1*jay] <= 0 | 931 is.na(wz[sind2, M1*jay]))) { 932 ind2[sind2[eim.kk.TF], jay] <- FALSE 933 } 934 935 936 lwr.ptr <- upr.ptr + 1 937 } # while 938 939 } # if 940 } # end of for (jay in 1:NOS) 941 942 943 944 945 946 947 948 for (jay in 1:NOS) { 949 run.varcov <- 0 950 ii.TF <- !ind2[, jay] # Not assigned above 951 if (any(ii.TF)) { 952 kkvec <- kmat[ii.TF, jay] 953 muvec <- munb[ii.TF, jay] 954 for (ii in 1:( .nsimEIM )) { 955 ysim <- rgaitnbinom(sum(ii.TF), kkvec, munb.p = muvec, 956 truncate = 0) 957 dl.dk <- digamma(ysim + kkvec) - digamma(kkvec) - 958 (ysim - muvec) / (muvec + kkvec) + 959 log1p(-muvec / (kkvec + muvec)) + 960 df0.dkmat[ii.TF, jay] / oneminusf0[ii.TF, jay] 961 run.varcov <- run.varcov + dl.dk^2 962 } # end of for loop 963 964 run.varcov <- c(run.varcov / .nsimEIM ) 965 ned2l.dk2 <- if (intercept.only) mean(run.varcov) else run.varcov 966 967 wz[ii.TF, M1*jay] <- ned2l.dk2 # * (dsize.deta[ii.TF, jay])^2 968 } 969 } # jay 970 971 972 973 wz[, M1*(1:NOS) ] <- wz[, M1*(1:NOS) ] * dsize.deta^2 974 975 976 977 978 979 980 save.weights <- !all(ind2) 981 982 ned2l.dmunb2 <- mymu / munb^2 - 983 ((1 + mymu/kmat) / kmat) / (1 + munb/kmat)^2 - 984 df02.dmunb2 / oneminusf0 - 985 (df0.dmunb / oneminusf0)^2 986 wz[, M1*(1:NOS) - 1] <- ned2l.dmunb2 * dmunb.deta^2 987 988 989 ned2l.dmunbsize <- (munb - mymu) / (munb + kmat)^2 - 990 df02.dkmat.dmunb / oneminusf0 - 991 df0.dmunb * df0.dkmat / oneminusf0^2 992 wz[, M + M1*(1:NOS) - 1] <- ned2l.dmunbsize * 993 dmunb.deta * dsize.deta 994 995 996 997 998 w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = NOS) 999 }), list( .cutoff.prob = cutoff.prob, .eps.trig = eps.trig, 1000 .max.support = max.support, 1001 .max.chunk.MB = max.chunk.MB, 1002 .nsimEIM = nsimEIM )))) 1003 1004} # posnegbinomial 1005 1006 1007 1008 1009 1010 1011dposgeom <- function(x, prob, log = FALSE) { 1012 dgeom(x - 1, prob = prob, log = log) 1013} 1014 1015 1016 1017pposgeom <- function(q, prob) { 1018 L <- max(length(q), length(prob)) 1019 if (length(q) != L) q <- rep_len(q, L) 1020 if (length(prob) != L) prob <- rep_len(prob, L) 1021 1022 ans <- ifelse(q < 1, 0, (pgeom(q, prob) - dgeom(0, prob)) 1023 / pgeom(0, prob, lower.tail = FALSE)) 1024 ans[prob == 1] <- NaN 1025 ans[prob == 0] <- NaN 1026 ans 1027 1028} 1029 1030 1031 1032qposgeom <- function(p, prob) { 1033 1034 1035 1036 1037 ans <- qgeom(pgeom(0, prob, lower.tail = FALSE) * p + dgeom(0, prob), 1038 prob) 1039 1040 ans[p == 1] <- Inf 1041 1042 ans[p <= 0] <- NaN 1043 ans[1 < p] <- NaN 1044 ans[prob == 0] <- NaN 1045 ans[prob == 1] <- NaN 1046 ans 1047} 1048 1049 1050 1051rposgeom <- function(n, prob) { 1052 ans <- qgeom(p = runif(n, min = dgeom(0, prob)), prob) 1053 ans[prob == 0] <- NaN 1054 ans[prob == 1] <- NaN 1055 ans 1056} 1057 1058 1059 1060 1061 1062 1063 pospoisson <- 1064 function(link = "loglink", 1065 type.fitted = c("mean", "lambda", "prob0"), 1066 expected = TRUE, 1067 ilambda = NULL, imethod = 1, zero = NULL, 1068 gt.1 = FALSE) { 1069 1070 link <- as.list(substitute(link)) 1071 earg <- link2list(link) 1072 link <- attr(earg, "function.name") 1073 1074 1075 if (!is.logical(expected) || length(expected) != 1) 1076 stop("bad input for argument 'expected'") 1077 if (length( ilambda) && !is.Numeric(ilambda, positive = TRUE)) 1078 stop("bad input for argument 'ilambda'") 1079 1080 type.fitted <- match.arg(type.fitted, 1081 c("mean", "lambda", "prob0"))[1] 1082 1083 1084 1085 1086 new("vglmff", 1087 blurb = c("Positive-Poisson distribution\n\n", 1088 "Links: ", 1089 namesof("lambda", link, earg = earg, tag = FALSE)), 1090 constraints = eval(substitute(expression({ 1091 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1092 predictors.names = predictors.names, 1093 M1 = 1) 1094 }), list( .zero = zero ))), 1095 1096 infos = eval(substitute(function(...) { 1097 list(M1 = 1, 1098 Q1 = 1, 1099 expected = TRUE, 1100 multipleResponses = TRUE, 1101 parameters.names = c("lambda"), 1102 link = .link , 1103 type.fitted = .type.fitted , 1104 expected = .expected , 1105 earg = .earg) 1106 }, list( .link = link, .earg = earg, 1107 .expected = expected, 1108 .type.fitted = type.fitted ))), 1109 1110 initialize = eval(substitute(expression({ 1111 temp5 <- 1112 w.y.check(w = w, y = y, 1113 Is.positive.y = TRUE, 1114 Is.integer.y = TRUE, 1115 ncol.w.max = Inf, 1116 ncol.y.max = Inf, 1117 out.wy = TRUE, 1118 colsyperw = 1, 1119 maximize = TRUE) 1120 w <- temp5$w 1121 y <- temp5$y 1122 1123 ncoly <- ncol(y) 1124 M1 <- 1 1125 extra$ncoly <- ncoly 1126 extra$M1 <- M1 1127 M <- M1 * ncoly 1128 extra$type.fitted <- .type.fitted 1129 extra$colnames.y <- colnames(y) 1130 1131 mynames1 <- param.names("lambda", ncoly, skip1 = TRUE) 1132 predictors.names <- namesof(mynames1, .link , earg = .earg, 1133 tag = FALSE) 1134 1135 if (!length(etastart)) { 1136 lambda.init <- if (length( .ilambda )) 1137 rep( .ilambda , length = n) else 1138 Init.mu(y = y, w = w, imethod = .imethod , 1139 imu = .ilambda ) 1140 etastart <- theta2eta(lambda.init, .link , earg = .earg ) 1141 } 1142 }), list( .link = link, .earg = earg, 1143 .ilambda = ilambda, .imethod = imethod, 1144 .type.fitted = type.fitted ))), 1145 linkinv = eval(substitute(function(eta, extra = NULL) { 1146 NOS <- NCOL(eta) / c(M1 = 1) 1147 type.fitted <- if (length(extra$type.fitted)) 1148 extra$type.fitted else { 1149 warning("cannot find 'type.fitted'. Returning the 'mean'.") 1150 "mean" 1151 } 1152 1153 type.fitted <- match.arg(type.fitted, 1154 c("mean", "lambda", "prob0"))[1] 1155 1156 lambda <- eta2theta(eta, .link , earg = .earg ) 1157 ans <- switch(type.fitted, 1158 "mean" = -lambda / expm1(-lambda), 1159 "lambda" = lambda, 1160 "prob0" = exp(-lambda)) # P(Y=0) as it were 1161 label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS) 1162 }, list( .link = link, .earg = earg ))), 1163 last = eval(substitute(expression({ 1164 misc$link <- rep_len( .link , M) 1165 names(misc$link) <- mynames1 1166 1167 misc$earg <- vector("list", M) 1168 names(misc$earg) <- mynames1 1169 for (ii in 1:M) 1170 misc$earg[[ii]] <- .earg 1171 1172 misc$M1 <- M1 1173 misc$expected <- TRUE 1174 misc$multipleResponses <- TRUE 1175 }), list( .link = link, .earg = earg, .expected = expected ))), 1176 loglikelihood = eval(substitute( 1177 function(mu, y, w, residuals = FALSE, eta, 1178 extra = NULL, summation = TRUE) { 1179 lambda <- eta2theta(eta, .link , earg = .earg ) 1180 if (residuals) { 1181 stop("loglikelihood residuals not implemented yet") 1182 } else { 1183 ll.elts <- c(w) * dgaitpois(y, lambda, truncate = 0, log = TRUE) 1184 if (summation) { 1185 sum(ll.elts) 1186 } else { 1187 ll.elts 1188 } 1189 } 1190 }, list( .link = link, .earg = earg ))), 1191 vfamily = c("pospoisson"), 1192 validparams = eval(substitute(function(eta, y, extra = NULL) { 1193 lambda <- eta2theta(eta, .link , earg = .earg ) 1194 lower.bound <- if ( .gt.1 ) 1 else 0 1195 okay1 <- all(is.finite(lambda)) && 1196 all(lower.bound < lambda) 1197 okay1 1198 }, list( .link = link, .earg = earg, .gt.1 = gt.1 ))), 1199 1200 1201 simslot = eval(substitute( 1202 function(object, nsim) { 1203 1204 pwts <- if (length(pwts <- object@prior.weights) > 0) 1205 pwts else weights(object, type = "prior") 1206 if (any(pwts != 1)) 1207 warning("ignoring prior weights") 1208 eta <- predict(object) 1209 lambda <- eta2theta(eta, .link , earg = .earg ) 1210 rgaitpois(nsim * length(lambda), lambda, truncate = 0) 1211 }, list( .link = link, .earg = earg ))), 1212 1213 1214 1215 1216 deriv = eval(substitute(expression({ 1217 lambda <- eta2theta(eta, .link , earg = .earg ) 1218 1219 temp6 <- expm1(lambda) 1220 dl.dlambda <- y / lambda - 1 - 1 / temp6 1221 1222 dlambda.deta <- dtheta.deta(lambda, .link , earg = .earg ) 1223 c(w) * dl.dlambda * dlambda.deta 1224 }), list( .link = link, .earg = earg ))), 1225 weight = eval(substitute(expression({ 1226 if ( .expected ) { 1227 ned2l.dlambda2 <- (1 + 1 / temp6) * (1/lambda - 1/temp6) 1228 wz <- ned2l.dlambda2 * dlambda.deta^2 1229 } else { 1230 d2l.dlambda2 <- y / lambda^2 - (1 + 1 / temp6 + 1) / temp6 1231 d2lambda.deta2 <- d2theta.deta2(lambda, .link , earg = .earg) 1232 wz <- (dlambda.deta^2) * d2l.dlambda2 - 1233 dl.dlambda * d2lambda.deta2 1234 } 1235 c(w) * wz 1236 }), list( .link = link, .earg = earg, .expected = expected )))) 1237} # pospoisson 1238 1239 1240 1241 1242 1243 1244 1245 1246 posbinomial <- 1247 function(link = "logitlink", 1248 multiple.responses = FALSE, parallel = FALSE, 1249 omit.constant = FALSE, 1250 1251 p.small = 1e-4, no.warning = FALSE, 1252 1253 zero = NULL) { 1254 1255 1256 1257 1258 link <- as.list(substitute(link)) 1259 earg <- link2list(link) 1260 link <- attr(earg, "function.name") 1261 1262 1263 1264 if (!is.logical(multiple.responses) || length(multiple.responses) != 1) 1265 stop("bad input for argument 'multiple.responses'") 1266 1267 if (!is.logical(omit.constant) || length(omit.constant) != 1) 1268 stop("bad input for argument 'omit.constant'") 1269 1270 1271 1272 if (!is.Numeric(p.small, positive = TRUE, length.arg = 1)) 1273 stop("bad input for argument 'p.small'") 1274 1275 1276 new("vglmff", 1277 blurb = c("Positive-binomial distribution\n\n", 1278 "Links: ", 1279 if (multiple.responses) 1280 c(namesof("prob1", link, earg = earg, tag = FALSE), 1281 ",...,", 1282 namesof("probM", link, earg = earg, tag = FALSE)) else 1283 namesof("prob", link, earg = earg, tag = FALSE), 1284 "\n"), 1285 constraints = eval(substitute(expression({ 1286 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1287 bool = .parallel , 1288 constraints = constraints) 1289 1290 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1291 predictors.names = predictors.names, 1292 M1 = 1) 1293 }), list( .parallel = parallel, .zero = zero ))), 1294 infos = eval(substitute(function(...) { 1295 list(M1 = 1, 1296 Q1 = 1, 1297 expected = TRUE, 1298 multipleResponses = .multiple.responses , 1299 parameters.names = c("prob"), 1300 p.small = .p.small , 1301 no.warning = .no.warning , 1302 zero = .zero ) 1303 }, list( .zero = zero, 1304 .p.small = p.small, 1305 .multiple.responses = multiple.responses, 1306 .no.warning = no.warning ))), 1307 1308 initialize = eval(substitute(expression({ 1309 1310 mustart.orig <- mustart 1311 if ( .multiple.responses ) { 1312 temp5 <- 1313 w.y.check(w = w, y = y, 1314 Is.positive.y = TRUE, 1315 ncol.w.max = Inf, 1316 ncol.y.max = Inf, 1317 out.wy = TRUE, 1318 colsyperw = 1, 1319 maximize = TRUE) 1320 w <- temp5$w 1321 y <- temp5$y 1322 1323 1324 ncoly <- ncol(y) 1325 M1 <- 1 1326 extra$ncoly <- ncoly 1327 extra$M1 <- M1 1328 M <- M1 * ncoly 1329 1330 1331 extra$p.small <- .p.small 1332 extra$no.warning <- .no.warning 1333 1334 extra$orig.w <- w 1335 mustart <- matrix(colSums(y) / colSums(w), # Not colSums(y * w)... 1336 n, ncoly, byrow = TRUE) 1337 1338 } else { 1339 eval(binomialff(link = .earg , # earg = .earg , 1340 earg.link = TRUE)@initialize) 1341 } 1342 1343 1344 if ( .multiple.responses ) { 1345 1346 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 1347 dn2 <- if (length(dn2)) { 1348 paste("E[", dn2, "]", sep = "") 1349 } else { 1350 param.names("prob", M) 1351 } 1352 predictors.names <- 1353 namesof(if (M > 1) dn2 else "prob", 1354 .link , earg = .earg, short = TRUE) 1355 1356 w <- matrix(w, n, ncoly) 1357 y <- y / w # Now sample proportion 1358 } else { 1359 predictors.names <- 1360 namesof("prob", .link , earg = .earg , tag = FALSE) 1361 } 1362 1363 if (length(extra)) extra$w <- w else extra <- list(w = w) 1364 1365 if (!length(etastart)) { 1366 mustart.use <- if (length(mustart.orig)) mustart.orig else mustart 1367 etastart <- cbind(theta2eta(mustart.use, .link , earg = .earg )) 1368 } 1369 mustart <- NULL 1370 1371 1372 1373 nvec <- if (NCOL(y) > 1) { 1374 NULL 1375 } else { 1376 if (is.numeric(extra$orig.w)) round(w / extra$orig.w) else 1377 round(w) 1378 } 1379 extra$tau <- if (length(nvec) && length(unique(nvec) == 1)) 1380 nvec[1] else NULL 1381 }), list( .link = link, 1382 .p.small = p.small, 1383 .no.warning = no.warning, 1384 .earg = earg, .multiple.responses = multiple.responses ))), 1385 1386 linkinv = eval(substitute(function(eta, extra = NULL) { 1387 w <- extra$w 1388 binprob <- eta2theta(eta, .link , earg = .earg ) 1389 nvec <- if ( .multiple.responses ) { 1390 w 1391 } else { 1392 if (is.numeric(extra$orig.w)) round(w / extra$orig.w) else 1393 round(w) 1394 } 1395 binprob / (1.0 - (1.0 - binprob)^nvec) 1396 }, 1397 1398 list( .link = link, .earg = earg, 1399 .multiple.responses = multiple.responses ))), 1400 last = eval(substitute(expression({ 1401 1402 1403 misc$link <- rep_len( .link , M) 1404 names(misc$link) <- if (M > 1) dn2 else "prob" 1405 1406 misc$earg <- vector("list", M) 1407 names(misc$earg) <- names(misc$link) 1408 for (ii in 1:M) 1409 misc$earg[[ii]] <- .earg 1410 1411 misc$expected <- TRUE 1412 misc$omit.constant <- .omit.constant 1413 misc$needto.omit.constant <- TRUE # Safety mechanism 1414 1415 1416 misc$multiple.responses <- .multiple.responses 1417 w <- as.numeric(w) 1418 1419 1420 1421 if (length(extra$tau)) { 1422 R <- tfit$qr$qr[1:ncol.X.vlm, 1:ncol.X.vlm, drop = FALSE] 1423 R[lower.tri(R)] <- 0 1424 tmp6 <- N.hat.posbernoulli(eta = eta, link = .link , 1425 earg = .earg , R = R, w = w, 1426 X.vlm = X.vlm.save, 1427 Hlist = Hlist, # 20150428; bug fixed here 1428 extra = extra, model.type = "0") 1429 extra$N.hat <- tmp6$N.hat 1430 extra$SE.N.hat <- tmp6$SE.N.hat 1431 } 1432 1433 1434 }), list( .link = link, .earg = earg, 1435 .multiple.responses = multiple.responses, 1436 .omit.constant = omit.constant ))), 1437 1438 loglikelihood = eval(substitute( 1439 function(mu, y, w, residuals = FALSE, eta, 1440 extra = NULL, 1441 summation = TRUE) { 1442 1443 ycounts <- if ( .multiple.responses ) { 1444 round(y * extra$orig.w) 1445 } else { 1446 if (is.numeric(extra$orig.w)) 1447 y * w / extra$orig.w else 1448 y * w # Convert proportions to counts 1449 } 1450 nvec <- if ( .multiple.responses ) { 1451 w 1452 } else { 1453 if (is.numeric(extra$orig.w)) 1454 round(w / extra$orig.w) else round(w) 1455 } 1456 use.orig.w <- if (is.numeric(extra$orig.w)) extra$orig.w else 1 1457 binprob <- eta2theta(eta, .link , earg = .earg ) 1458 1459 if (residuals) { 1460 stop("loglikelihood residuals not implemented yet") 1461 } else { 1462 answer <- c(use.orig.w) * 1463 dgaitbinom(ycounts, nvec, binprob, truncate = 0, log = TRUE) 1464 if ( .omit.constant ) { 1465 answer <- answer - c(use.orig.w) * lchoose(nvec, ycounts) 1466 } 1467 ll.elts <- answer 1468 if (summation) { 1469 sum(ll.elts) 1470 } else { 1471 ll.elts 1472 } 1473 } 1474 }, list( .link = link, .earg = earg, 1475 .multiple.responses = multiple.responses, 1476 .omit.constant = omit.constant ))), 1477 1478 vfamily = c("posbinomial"), 1479 validparams = eval(substitute(function(eta, y, extra = NULL) { 1480 binprob <- eta2theta(eta, .link , earg = .earg ) 1481 okay1 <- all(is.finite(binprob)) && all(0 < binprob & binprob < 1) 1482 okay1 1483 }, list( .link = link, .earg = earg ))), 1484 1485 1486 1487 1488 1489 simslot = eval(substitute( 1490 function(object, nsim) { 1491 1492 pwts <- if (length(pwts <- object@prior.weights) > 0) 1493 pwts else weights(object, type = "prior") 1494 1495 if ( .multiple.responses ) 1496 stop("cannot run simulate() when 'multiple.responses = TRUE'") 1497 1498 eta <- predict(object) 1499 binprob <- eta2theta(eta, .link , earg = .earg ) 1500 1501 extra <- object@extra 1502 w <- extra$w # Usual code 1503 w <- pwts # 20140101 1504 1505 1506 nvec <- if ( .multiple.responses ) { 1507 w 1508 } else { 1509 if (is.numeric(extra$orig.w)) round(w / extra$orig.w) else 1510 round(w) 1511 } 1512 rgaitbinom(nsim * length(eta), nvec, binprob, truncate = 0) 1513 }, list( .link = link, .earg = earg, 1514 .multiple.responses = multiple.responses, 1515 .omit.constant = omit.constant ))), 1516 1517 1518 1519 1520 1521 deriv = eval(substitute(expression({ 1522 use.orig.w <- if (is.numeric(extra$orig.w)) extra$orig.w else 1523 rep_len(1, n) 1524 1525 nvec <- if ( .multiple.responses ) { 1526 w 1527 } else { 1528 if (is.numeric(extra$orig.w)) round(w / extra$orig.w) else 1529 round(w) 1530 } 1531 binprob <- eta2theta(eta, .link , earg = .earg ) 1532 dmu.deta <- dtheta.deta(binprob, .link , earg = .earg ) 1533 1534 temp1 <- 1 - (1 - binprob)^nvec 1535 temp2 <- (1 - binprob)^2 1536 temp3 <- (1 - binprob)^(nvec-2) 1537 1538 dl.dmu <- y / binprob - (1 - y) / (1 - binprob) - 1539 (1 - binprob) * temp3 / temp1 1540 1541 c(w) * dl.dmu * dmu.deta 1542 }), list( .link = link, .earg = earg, 1543 .multiple.responses = multiple.responses ))), 1544 weight = eval(substitute(expression({ 1545 1546 ned2l.dmu2 <- 1 / (binprob * temp1) + 1547 (1 - mu) / temp2 - 1548 (nvec-1) * temp3 / temp1 - 1549 nvec * (temp2^(nvec-1)) / temp1^2 1550 1551 1552 1553 wz <- c(w) * ned2l.dmu2 * dmu.deta^2 1554 wz 1555 }), list( .link = link, .earg = earg, 1556 .multiple.responses = multiple.responses )))) 1557} # posbinomial 1558 1559 1560 1561 1562 1563 1564 posbernoulli.t <- 1565 function(link = "logitlink", 1566 1567 parallel.t = FALSE ~ 1, 1568 1569 1570 1571 iprob = NULL, 1572 1573 p.small = 1e-4, no.warning = FALSE, 1574 type.fitted = c("probs", "onempall0")) { 1575 1576 1577 1578 1579 1580 1581 type.fitted <- match.arg(type.fitted, 1582 c("probs", "onempall0"))[1] 1583 1584 1585 apply.parint <- FALSE 1586 1587 1588 1589 link <- as.list(substitute(link)) 1590 earg <- link2list(link) 1591 link <- attr(earg, "function.name") 1592 1593 1594 if (length(iprob)) 1595 if (!is.Numeric(iprob, positive = TRUE) || 1596 max(iprob) >= 1) 1597 stop("argument 'iprob' must have values in (0, 1)") 1598 1599 if (!is.logical(apply.parint) || 1600 length(apply.parint) != 1) 1601 stop("argument 'apply.parint' must be a single logical") 1602 1603 if (!is.Numeric(p.small, positive = TRUE, length.arg = 1)) 1604 stop("bad input for argument 'p.small'") 1605 1606 1607 new("vglmff", 1608 blurb = c("Positive-Bernoulli (capture-recapture) model ", 1609 "with temporal effects (M_{t}/M_{th})\n\n", 1610 "Links: ", 1611 namesof("prob1", link, earg = earg, tag = FALSE), ", ", 1612 namesof("prob2", link, earg = earg, tag = FALSE), ", ..., ", 1613 namesof("probM", link, earg = earg, tag = FALSE), 1614 "\n"), 1615 constraints = eval(substitute(expression({ 1616 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1617 bool = .parallel.t , 1618 constraints = constraints, 1619 apply.int = .apply.parint , # TRUE, 1620 cm.default = diag(M), 1621 cm.intercept.default = diag(M)) 1622 }), list( .parallel.t = parallel.t, 1623 .apply.parint = apply.parint ))), 1624 infos = eval(substitute(function(...) { 1625 list(M1 = 1, 1626 Q1 = NA, 1627 expected = TRUE, 1628 multipleResponses = TRUE, 1629 parameters.names = c("prob"), 1630 p.small = .p.small , 1631 type.fitted = .type.fitted , 1632 no.warning = .no.warning , 1633 apply.parint = .apply.parint , 1634 parallel.t = .parallel.t ) 1635 }, list( .parallel.t = parallel.t, 1636 .p.small = p.small, 1637 .no.warning = no.warning, 1638 .type.fitted = type.fitted, 1639 .apply.parint = apply.parint ))), 1640 1641 initialize = eval(substitute(expression({ 1642 M1 <- 1 1643 1644 mustart.orig <- mustart 1645 y <- as.matrix(y) 1646 M <- ncoly <- ncol(y) 1647 extra$ncoly <- ncoly <- ncol(y) 1648 extra$tau <- tau <- ncol(y) 1649 extra$orig.w <- w 1650 1651 extra$p.small <- .p.small 1652 extra$no.warning <- .no.warning 1653 1654 extra$type.fitted <- .type.fitted 1655 extra$colnames.y <- colnames(y) 1656 1657 w <- matrix(w, n, ncoly) 1658 mustart <- matrix(colSums(y) / colSums(w), 1659 n, ncol(y), byrow = TRUE) 1660 mustart[mustart == 0] <- 0.05 1661 mustart[mustart == 1] <- 0.95 1662 1663 if (ncoly == 1) 1664 stop("the response is univariate, therefore use posbinomial()") 1665 1666 1667 1668 1669 1670 1671 if (!all(y == 0 | y == 1)) 1672 stop("response must contain 0s and 1s only") 1673 if (!all(w == 1)) 1674 stop("argument 'weight' must contain 1s only") 1675 1676 1677 1678 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 1679 dn2 <- if (length(dn2)) { 1680 paste("E[", dn2, "]", sep = "") 1681 } else { 1682 param.names("prob", M) 1683 } 1684 1685 1686 predictors.names <- namesof(dn2, .link , earg = .earg, short = TRUE) 1687 1688 1689 if (length(extra)) extra$w <- w else extra <- list(w = w) 1690 1691 if (!length(etastart)) { 1692 mustart.use <- if (length(mustart.orig)) { 1693 mustart.orig 1694 } else { 1695 mustart 1696 } 1697 etastart <- cbind(theta2eta(mustart.use, .link , earg = .earg )) 1698 } 1699 mustart <- NULL 1700 }), list( .link = link, .earg = earg, 1701 .p.small = p.small, 1702 .type.fitted = type.fitted, 1703 .no.warning = no.warning 1704 ))), 1705 linkinv = eval(substitute(function(eta, extra = NULL) { 1706 type.fitted <- 1707 if (length(extra$type.fitted)) extra$type.fitted else { 1708 warning("cannot find 'type.fitted'. Returning the 'probs'.") 1709 "probs" 1710 } 1711 1712 type.fitted <- match.arg(type.fitted, 1713 c("probs", "onempall0"))[1] 1714 1715 tau <- extra$ncoly 1716 probs <- eta2theta(eta, .link , earg = .earg ) 1717 logAA0 <- rowSums(log1p(-probs)) 1718 AA0 <- exp(logAA0) 1719 AAA <- exp(log1p(-AA0)) # 1 - AA0 1720 1721 1722 1723 fv <- probs / AAA 1724 ans <- switch(type.fitted, 1725 "probs" = fv, 1726 "onempall0" = AAA) 1727 label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS) 1728}, list( .link = link, .earg = earg ))), 1729 last = eval(substitute(expression({ 1730 extra$w <- NULL # Kill it off 1731 1732 1733 misc$link <- rep_len( .link , M) 1734 names(misc$link) <- if (M > 1) dn2 else "prob" 1735 1736 misc$earg <- vector("list", M) 1737 names(misc$earg) <- names(misc$link) 1738 for (ii in 1:M) misc$earg[[ii]] <- .earg 1739 1740 1741 misc$multiple.responses <- TRUE 1742 misc$iprob <- .iprob 1743 1744 1745 R <- tfit$qr$qr[1:ncol.X.vlm, 1:ncol.X.vlm, drop = FALSE] 1746 R[lower.tri(R)] <- 0 1747 tmp6 <- N.hat.posbernoulli(eta = eta, link = .link , earg = .earg , 1748 R = R, w = w, 1749 X.vlm = X.vlm.save, 1750 Hlist = Hlist, # 20150428 bug fixed here 1751 extra = extra, model.type = "t") 1752 extra$N.hat <- tmp6$N.hat 1753 extra$SE.N.hat <- tmp6$SE.N.hat 1754 1755 1756 1757 1758 misc$parallel.t <- .parallel.t 1759 misc$apply.parint <- .apply.parint 1760 }), list( .link = link, .earg = earg, 1761 .parallel.t = parallel.t, 1762 .apply.parint = apply.parint, 1763 .iprob = iprob ))), 1764 loglikelihood = eval(substitute( 1765 function(mu, y, w, residuals = FALSE, eta, 1766 extra = NULL, 1767 summation = TRUE) { 1768 1769 ycounts <- y 1770 use.orig.w <- if (length(extra$orig.w)) extra$orig.w else 1 1771 1772 probs <- eta2theta(eta, .link , earg = .earg ) 1773 1774 if (residuals) { 1775 stop("loglikelihood residuals not implemented yet") 1776 } else { 1777 1778 1779 ll.elts <- 1780 c(use.orig.w) * 1781 dposbern(x = ycounts, # size = 1, # Bernoulli trials 1782 prob = probs, prob0 = probs, log = TRUE) 1783 if (summation) { 1784 sum(ll.elts) 1785 } else { 1786 ll.elts 1787 } 1788 } 1789 }, list( .link = link, .earg = earg ))), 1790 vfamily = c("posbernoulli.t"), 1791 validparams = eval(substitute(function(eta, y, extra = NULL) { 1792 probs <- eta2theta(eta, .link , earg = .earg ) 1793 okay1 <- all(is.finite(probs)) && all(0 < probs & probs < 1) 1794 okay1 1795 }, list( .link = link, .earg = earg ))), 1796 1797 1798 deriv = eval(substitute(expression({ 1799 probs <- eta2theta(eta, .link , earg = .earg ) 1800 dprobs.deta <- dtheta.deta(probs, .link , earg = .earg ) 1801 1802 logAA0 <- rowSums(log1p(-probs)) 1803 AA0 <- exp(logAA0) 1804 AAA <- exp(log1p(-AA0)) # 1 - AA0 1805 1806 B.s <- AA0 / (1 - probs) 1807 B.st <- array(AA0, c(n, M, M)) 1808 for (slocal in 1:(M-1)) 1809 for (tlocal in (slocal+1):M) 1810 B.st[, slocal, tlocal] <- 1811 B.st[, tlocal, slocal] <- B.s[, slocal] / (1 - probs[, tlocal]) 1812 1813 temp2 <- (1 - probs)^2 1814 dl.dprobs <- y / probs - (1 - y) / (1 - probs) - B.s / AAA 1815 1816 deriv.ans <- c(w) * dl.dprobs * dprobs.deta 1817 deriv.ans 1818 }), list( .link = link, .earg = earg ))), 1819 weight = eval(substitute(expression({ 1820 1821 ned2l.dprobs2 <- 1 / (probs * AAA) + 1 / temp2 - 1822 probs / (AAA * temp2) - (B.s / AAA)^2 1823 1824 wz <- matrix(NA_real_, n, dimm(M)) 1825 wz[, 1:M] <- ned2l.dprobs2 * (dprobs.deta^2) 1826 1827 for (slocal in 1:(M-1)) 1828 for (tlocal in (slocal+1):M) 1829 wz[, iam(slocal, tlocal, M = M)] <- 1830 dprobs.deta[, slocal] * dprobs.deta[, tlocal] * 1831 (B.st[, slocal, tlocal] + 1832 B.s [, slocal] * B.s [, tlocal] / AAA) / (-AAA) 1833 wz 1834 }), list( .link = link, .earg = earg )))) 1835} # posbernoulli.t 1836 1837 1838 1839 1840 1841 1842 posbernoulli.b <- 1843 function(link = "logitlink", 1844 1845 1846 drop.b = FALSE ~ 1, 1847 1848 1849 type.fitted = c("likelihood.cond", "mean.uncond"), 1850 1851 I2 = FALSE, 1852 ipcapture = NULL, 1853 iprecapture = NULL, 1854 p.small = 1e-4, no.warning = FALSE) { 1855 1856 1857 1858 1859 type.fitted <- match.arg(type.fitted, 1860 c("likelihood.cond", "mean.uncond"))[1] 1861 1862 link <- as.list(substitute(link)) 1863 earg <- link2list(link) 1864 link <- attr(earg, "function.name") 1865 1866 apply.parint.b <- FALSE 1867 1868 1869 if (length(ipcapture)) 1870 if (!is.Numeric(ipcapture, positive = TRUE) || 1871 max(ipcapture) >= 1) 1872 stop("argument 'ipcapture' must have values in (0, 1)") 1873 if (length(iprecapture)) 1874 if (!is.Numeric(iprecapture, positive = TRUE) || 1875 max(iprecapture) >= 1) 1876 stop("argument 'iprecapture' must have values in (0, 1)") 1877 1878 if (!is.logical(I2) || 1879 length(I2) != 1) 1880 stop("argument 'I2' must be a single logical") 1881 1882 1883 if (!is.Numeric(p.small, positive = TRUE, length.arg = 1)) 1884 stop("bad input for argument 'p.small'") 1885 1886 1887 1888 1889 new("vglmff", 1890 blurb = c("Positive-Bernoulli (capture-recapture) model ", 1891 "with behavioural effects (M_{b}/M_{bh})\n\n", 1892 "Links: ", 1893 namesof("pcapture", link, earg = earg, tag = FALSE), ", ", 1894 namesof("precapture", link, earg = earg, tag = FALSE), 1895 "\n"), 1896 1897 constraints = eval(substitute(expression({ 1898 1899 cm.intercept.default <- if ( .I2 ) diag(2) else cbind(0:1, 1) 1900 1901 constraints <- cm.VGAM(matrix(1, 2, 1), x = x, 1902 bool = .drop.b , 1903 constraints = constraints, 1904 apply.int = .apply.parint.b , # TRUE, 1905 cm.default = cm.intercept.default, # diag(2), 1906 cm.intercept.default = cm.intercept.default) 1907 }), list( .drop.b = drop.b, 1908 .I2 = I2, 1909 .apply.parint.b = apply.parint.b ))), 1910 1911 infos = eval(substitute(function(...) { 1912 list(M1 = 2, 1913 expected = TRUE, 1914 multipleResponses = FALSE, 1915 parameters.names = c("pcapture", "precapture"), 1916 p.small = .p.small , 1917 no.warning = .no.warning , 1918 type.fitted = .type.fitted , 1919 apply.parint.b = .apply.parint.b ) 1920 }, list( 1921 .apply.parint.b = apply.parint.b, 1922 .p.small = p.small, 1923 .no.warning = no.warning, 1924 .type.fitted = type.fitted 1925 ))), 1926 1927 initialize = eval(substitute(expression({ 1928 M1 <- 2 1929 if (!is.matrix(y) || ncol(y) == 1) 1930 stop("the response appears to be univariate") 1931 1932 if (!all(y == 0 | y == 1)) 1933 stop("response must contain 0s and 1s only") 1934 1935 orig.y <- y 1936 extra$orig.w <- w 1937 extra$tau <- tau <- ncol(y) 1938 extra$ncoly <- ncoly <- ncol(y) 1939 extra$type.fitted <- .type.fitted 1940 extra$colnames.y <- colnames(y) 1941 1942 1943 extra$p.small <- .p.small 1944 extra$no.warning <- .no.warning 1945 1946 1947 1948 1949 mustart.orig <- mustart 1950 M <- 2 1951 1952 1953 tmp3 <- aux.posbernoulli.t(y, rename = FALSE) 1954 y0i <- extra$y0i <- tmp3$y0i 1955 yr0i <- extra$yr0i <- tmp3$yr0i 1956 yr1i <- extra$yr1i <- tmp3$yr1i 1957 cap1 <- extra$cap1 <- tmp3$cap1 1958 cap.hist1 <- extra$cap.hist1 <- tmp3$cap.hist1 1959 1960 1961 temp5 <- 1962 w.y.check(w = w, y = y, 1963 Is.integer.y = TRUE, 1964 Is.nonnegative.y = TRUE, 1965 ncol.w.max = 1, 1966 ncol.y.min = 2, 1967 ncol.y.max = Inf, 1968 out.wy = TRUE, 1969 colsyperw = ncol(y), 1970 maximize = TRUE) 1971 w <- temp5$w # Retain the 0-1 response 1972 y <- temp5$y # Retain the 0-1 response 1973 1974 mustart <- matrix(colMeans(y), n, tau, byrow = TRUE) 1975 mustart <- (mustart + orig.y) / 2 1976 1977 1978 1979 1980 predictors.names <- 1981 c(namesof( "pcapture", .link , earg = .earg, short = TRUE), 1982 namesof("precapture", .link , earg = .earg, short = TRUE)) 1983 1984 if (!length(etastart)) { 1985 mustart.use <- if (length(mustart.orig)) { 1986 mustart.orig 1987 } else { 1988 mustart 1989 } 1990 1991 etastart <- 1992 cbind(theta2eta(rowMeans(mustart.use), .link , earg = .earg ), 1993 theta2eta(rowMeans(mustart.use), .link , earg = .earg )) 1994 1995 if (length( .ipcapture )) 1996 etastart[, 1] <- theta2eta( .ipcapture , .link , earg = .earg ) 1997 if (length( .iprecapture )) 1998 etastart[, 2] <- theta2eta( .iprecapture , .link , earg = .earg ) 1999 } 2000 mustart <- NULL 2001 }), list( .link = link, .earg = earg, 2002 .type.fitted = type.fitted, 2003 .p.small = p.small, 2004 .no.warning = no.warning, 2005 .ipcapture = ipcapture, 2006 .iprecapture = iprecapture 2007 ))), 2008 linkinv = eval(substitute(function(eta, extra = NULL) { 2009 cap.probs <- eta2theta(eta[, 1], .link , earg = .earg ) 2010 rec.probs <- eta2theta(eta[, 2], .link , earg = .earg ) 2011 tau <- extra$tau 2012 prc <- matrix(cap.probs, nrow(eta), tau) 2013 prr <- matrix(rec.probs, nrow(eta), tau) 2014 logQQQ <- rowSums(log1p(-prc)) 2015 QQQ <- exp(logQQQ) 2016 AAA <- exp(log1p(-QQQ)) # 1 - QQQ 2017 2018 2019 type.fitted <- if (length(extra$type.fitted)) extra$type.fitted else { 2020 warning("cannot find 'type.fitted'. ", 2021 "Returning 'likelihood.cond'.") 2022 "likelihood.cond" 2023 } 2024 2025 2026 type.fitted <- match.arg(type.fitted, 2027 c("likelihood.cond", "mean.uncond"))[1] 2028 2029 2030 if ( type.fitted == "likelihood.cond") { 2031 probs.numer <- prr 2032 mat.index <- cbind(1:nrow(prc), extra$cap1) 2033 probs.numer[mat.index] <- prc[mat.index] 2034 probs.numer[extra$cap.hist1 == 0] <- prc[extra$cap.hist1 == 0] 2035 fv <- probs.numer / AAA 2036 } else { 2037 2038 fv <- prc - prr 2039 for (jay in 2:tau) 2040 fv[, jay] <- fv[, jay-1] * (1 - cap.probs) 2041 fv <- (fv + prr) / AAA 2042 } 2043 2044 label.cols.y(fv, colnames.y = extra$colnames.y, NOS = tau) 2045 }, list( .link = link, .earg = earg ))), 2046 last = eval(substitute(expression({ 2047 2048 misc$link <- c( .link , .link ) 2049 names(misc$link) <- predictors.names 2050 2051 misc$earg <- vector("list", M) 2052 names(misc$earg) <- names(misc$link) 2053 misc$earg[[1]] <- .earg 2054 misc$earg[[2]] <- .earg 2055 2056 misc$expected <- TRUE 2057 misc$multiple.responses <- TRUE 2058 misc$ipcapture <- .ipcapture 2059 misc$iprecapture <- .iprecapture 2060 misc$drop.b <- .drop.b 2061 misc$multipleResponses <- FALSE 2062 misc$apply.parint.b <- .apply.parint.b 2063 2064 2065 2066 R <- tfit$qr$qr[1:ncol.X.vlm, 1:ncol.X.vlm, drop = FALSE] 2067 R[lower.tri(R)] <- 0 2068 tmp6 <- N.hat.posbernoulli(eta = eta, link = .link , earg = .earg , 2069 R = R, w = w, 2070 X.vlm = X.vlm.save, 2071 Hlist = Hlist, # 20150428; bug fixed here 2072 extra = extra, model.type = "b") 2073 extra$N.hat <- tmp6$N.hat 2074 extra$SE.N.hat <- tmp6$SE.N.hat 2075 2076 2077 }), list( .link = link, .earg = earg, 2078 .drop.b = drop.b, 2079 .ipcapture = ipcapture, 2080 .iprecapture = iprecapture, 2081 .apply.parint.b = apply.parint.b 2082 ))), 2083 loglikelihood = eval(substitute( 2084 function(mu, y, w, residuals = FALSE, eta, 2085 extra = NULL, 2086 summation = TRUE) { 2087 2088 tau <- extra$ncoly 2089 ycounts <- y 2090 use.orig.w <- if (length(extra$orig.w)) extra$orig.w else 1 2091 2092 cap.probs <- eta2theta(eta[, 1], .link , earg = .earg ) 2093 rec.probs <- eta2theta(eta[, 2], .link , earg = .earg ) 2094 prc <- matrix(cap.probs, nrow(eta), tau) 2095 prr <- matrix(rec.probs, nrow(eta), tau) 2096 2097 if (residuals) { 2098 stop("loglikelihood residuals not implemented yet") 2099 } else { 2100 probs.numer <- prr 2101 mat.index <- cbind(1:nrow(prc), extra$cap1) 2102 probs.numer[mat.index] <- prc[mat.index] 2103 probs.numer[extra$cap.hist1 == 0] <- prc[extra$cap.hist1 == 0] 2104 2105 ll.elts <- 2106 c(use.orig.w) * 2107 dposbern(x = ycounts, # Bernoulli trials 2108 prob = probs.numer, prob0 = prc, log = TRUE) 2109 if (summation) { 2110 sum(ll.elts) 2111 } else { 2112 ll.elts 2113 } 2114 } 2115 }, list( .link = link, .earg = earg ))), 2116 vfamily = c("posbernoulli.b"), 2117 validparams = eval(substitute(function(eta, y, extra = NULL) { 2118 cap.probs <- eta2theta(eta[, 1], .link , earg = .earg ) 2119 rec.probs <- eta2theta(eta[, 2], .link , earg = .earg ) 2120 okay1 <- all(is.finite(cap.probs)) && 2121 all(0 < cap.probs & cap.probs < 1) && 2122 all(is.finite(rec.probs)) && 2123 all(0 < rec.probs & rec.probs < 1) 2124 okay1 2125 }, list( .link = link, .earg = earg ))), 2126 2127 2128 2129 2130 2131 2132 2133 2134 2135 2136 2137 2138 2139 deriv = eval(substitute(expression({ 2140 cap.probs <- eta2theta(eta[, 1], .link , earg = .earg ) 2141 rec.probs <- eta2theta(eta[, 2], .link , earg = .earg ) 2142 y0i <- extra$y0i 2143 yr0i <- extra$yr0i 2144 yr1i <- extra$yr1i 2145 cap1 <- extra$cap1 2146 tau <- extra$tau 2147 2148 dcapprobs.deta <- dtheta.deta(cap.probs, .link , earg = .earg ) 2149 drecprobs.deta <- dtheta.deta(rec.probs, .link , earg = .earg ) 2150 2151 QQQ <- (1 - cap.probs)^tau 2152 dl.dcap <- 1 / cap.probs - 2153 y0i / (1 - cap.probs) - 2154 tau * ((1 - cap.probs)^(tau - 1)) / (1 - QQQ) 2155 2156 dl.drec <- yr1i / rec.probs - 2157 yr0i / (1 - rec.probs) 2158 2159 2160 deriv.ans <- c(w) * cbind(dl.dcap * dcapprobs.deta, 2161 dl.drec * drecprobs.deta) 2162 deriv.ans 2163 }), list( .link = link, .earg = earg ))), 2164 2165 weight = eval(substitute(expression({ 2166 2167 wz <- matrix(0, n, M) # Diagonal EIM 2168 2169 2170 2171 dA.dcapprobs <- -tau * ((1 - QQQ) * (tau-1) * (1 - cap.probs)^(tau-2) + 2172 tau * (1 - cap.probs)^(2*tau -2)) / (1 - QQQ)^2 2173 2174 2175 2176 2177 2178 prc <- matrix(cap.probs, n, tau) 2179 prr <- matrix(rec.probs, n, tau) 2180 2181 dQ.dprc <- -QQQ / (1 - prc) 2182 QQQcummat <- exp(t( apply(log1p(-prc), 1, cumsum))) 2183 2184 2185 2186 GGG <- (1 - QQQ - cap.probs * (1 + (tau-1) * QQQ)) / ( 2187 cap.probs * (1-cap.probs)^2) 2188 wz.pc <- GGG / (1 - QQQ) + 1 / cap.probs^2 + dA.dcapprobs 2189 wz[, iam(1, 1, M = M)] <- wz.pc * dcapprobs.deta^2 # Efficient 2190 2191 2192 2193 2194 2195 wz.pr <- (tau - (1 - QQQ) / cap.probs) / ( 2196 rec.probs * (1 - rec.probs) * (1 - QQQ)) 2197 wz[, iam(2, 2, M = M)] <- wz.pr * drecprobs.deta^2 2198 2199 2200 2201 2202 wz <- c(w) * wz 2203 wz 2204 }), list( .link = link, .earg = earg )))) 2205} # posbernoulli.b 2206 2207 2208 2209 2210 2211 2212 posbernoulli.tb <- 2213 function(link = "logitlink", 2214 parallel.t = FALSE ~ 1, 2215 parallel.b = FALSE ~ 0, 2216 drop.b = FALSE ~ 1, 2217 type.fitted = c("likelihood.cond", "mean.uncond"), 2218 imethod = 1, 2219 iprob = NULL, 2220 p.small = 1e-4, no.warning = FALSE, 2221 ridge.constant = 0.0001, # 20181020 2222 ridge.power = -4) { 2223 2224 2225 2226 apply.parint.t <- FALSE 2227 apply.parint.b <- TRUE 2228 apply.parint.d <- FALSE # For 'drop.b' actually. 2229 2230 link <- as.list(substitute(link)) 2231 earg <- link2list(link) 2232 link <- attr(earg, "function.name") 2233 2234 type.fitted <- match.arg(type.fitted, 2235 c("likelihood.cond", "mean.uncond"))[1] 2236 2237 2238 if (!is.Numeric(imethod, length.arg = 1, 2239 integer.valued = TRUE, positive = TRUE) || 2240 imethod > 2) 2241 stop("argument 'imethod' must be 1 or 2") 2242 2243 2244 if (!is.Numeric(ridge.constant) || 2245 ridge.constant < 0) 2246 warning("argument 'ridge.constant' should be non-negative") 2247 if (!is.Numeric(ridge.power) || 2248 ridge.power > 0) 2249 warning("argument 'ridge.power' should be non-positive") 2250 2251 2252 if (length(iprob)) 2253 if (!is.Numeric(iprob, positive = TRUE) || 2254 max(iprob) >= 1) 2255 stop("argument 'iprob' must have values in (0, 1)") 2256 2257 2258 if (!is.Numeric(p.small, positive = TRUE, length.arg = 1)) 2259 stop("bad input for argument 'p.small'") 2260 2261 2262 2263 new("vglmff", 2264 blurb = c("Positive-Bernoulli (capture-recapture) model\n", 2265 "with temporal and behavioural effects (M_{tb}/M_{tbh})\n\n", 2266 "Links: ", 2267 namesof("pcapture.1", link, earg = earg, tag = FALSE), 2268 ", ..., ", 2269 namesof("pcapture.tau", link, earg = earg, tag = FALSE), 2270 ", ", 2271 namesof("precapture.2", link, earg = earg, tag = FALSE), 2272 ", ..., ", 2273 namesof("precapture.tau", link, earg = earg, tag = FALSE)), 2274 constraints = eval(substitute(expression({ 2275 2276 2277 constraints.orig <- constraints 2278 cm1.d <- 2279 cmk.d <- matrix(0, M, 1) # All 0s inside 2280 con.d <- cm.VGAM(matrix(1, M, 1), x = x, 2281 bool = .drop.b , 2282 constraints = constraints.orig, 2283 apply.int = .apply.parint.d , # FALSE, 2284 cm.default = cmk.d, 2285 cm.intercept.default = cm1.d) 2286 2287 2288 2289 cm1.t <- 2290 cmk.t <- rbind(diag(tau), diag(tau)[-1, ]) # More readable 2291 con.t <- cm.VGAM(matrix(1, M, 1), x = x, 2292 bool = .parallel.t , # Same as .parallel.b 2293 constraints = constraints.orig, 2294 apply.int = .apply.parint.t , # FALSE, 2295 cm.default = cmk.t, 2296 cm.intercept.default = cm1.t) 2297 2298 2299 2300 cm1.b <- 2301 cmk.b <- rbind(matrix(0, tau, tau-1), diag(tau-1)) 2302 con.b <- cm.VGAM(matrix(c(rep_len(0, tau ), 2303 rep_len(1, tau-1)), M, 1), x = x, 2304 bool = .parallel.b , # Same as .parallel.b 2305 constraints = constraints.orig, 2306 apply.int = .apply.parint.b , # FALSE, 2307 cm.default = cmk.b, 2308 cm.intercept.default = cm1.b) 2309 2310 con.use <- con.b 2311 for (klocal in seq_along(con.b)) { 2312 con.use[[klocal]] <- 2313 cbind(if (any(con.d[[klocal]] == 1)) NULL else con.b[[klocal]], 2314 con.t[[klocal]]) 2315 2316 } 2317 2318 2319 constraints <- con.use 2320 2321 }), list( .parallel.t = parallel.t, 2322 .parallel.b = parallel.b, 2323 .drop.b = drop.b, 2324 .apply.parint.b = apply.parint.b, 2325 .apply.parint.d = apply.parint.d, 2326 .apply.parint.t = apply.parint.t ))), 2327 infos = eval(substitute(function(...) { 2328 list(M1 = 2, 2329 expected = TRUE, 2330 multipleResponses = TRUE, 2331 parameters.names = as.character(NA), 2332 ridge.constant = .ridge.constant , 2333 ridge.power = .ridge.power , 2334 drop.b = .drop.b, 2335 imethod = .imethod , 2336 type.fitted = .type.fitted , 2337 p.small = .p.small , 2338 no.warning = .no.warning , 2339 apply.parint.b = .apply.parint.b , 2340 apply.parint.t = .apply.parint.t , 2341 parallel.t = .parallel.t , 2342 parallel.b = .parallel.b ) 2343 }, list( .parallel.t = parallel.t, 2344 .parallel.b = parallel.b, 2345 .drop.b = drop.b, 2346 .type.fitted = type.fitted, 2347 .p.small = p.small, 2348 .no.warning = no.warning, 2349 .imethod = imethod, 2350 .ridge.constant = ridge.constant, 2351 .ridge.power = ridge.power, 2352 .apply.parint.b = apply.parint.b, 2353 .apply.parint.t = apply.parint.t ))), 2354 2355 initialize = eval(substitute(expression({ 2356 M1 <- 2 # Not quite true 2357 2358 2359 if (NCOL(w) > 1) 2360 stop("variable 'w' should be a vector or one-column matrix") 2361 w <- c(w) # Make it a vector 2362 2363 mustart.orig <- mustart 2364 y <- as.matrix(y) 2365 extra$tau <- tau <- ncol(y) 2366 extra$ncoly <- ncoly <- ncol(y) 2367 extra$orig.w <- w 2368 extra$ycounts <- y 2369 extra$type.fitted <- .type.fitted 2370 extra$colnames.y <- colnames(y) 2371 M <- M1 * tau - 1 # recap.prob.1 is unused 2372 2373 2374 mustart <- (y + matrix(apply(y, 2, weighted.mean, w = w), 2375 n, tau, byrow = TRUE)) / 2 2376 mustart[mustart < 0.01] <- 0.01 2377 mustart[mustart > 0.99] <- 0.99 2378 2379 mustart <- cbind(mustart, mustart[, -1]) 2380 2381 2382 2383 2384 extra$p.small <- .p.small 2385 extra$no.warning <- .no.warning 2386 2387 2388 2389 2390 2391 if (!all(y == 0 | y == 1)) 2392 stop("response must contain 0s and 1s only") 2393 2394 2395 tmp3 <- aux.posbernoulli.t(y) 2396 cap.hist1 <- extra$cap.hist1 <- tmp3$cap.hist1 2397 2398 2399 dn2.cap <- param.names("pcapture.", ncoly) 2400 dn2.recap <- param.names("precapture.", ncoly)[-1] 2401 2402 predictors.names <- c( 2403 namesof(dn2.cap, .link , earg = .earg, short = TRUE), 2404 namesof(dn2.recap, .link , earg = .earg, short = TRUE)) 2405 2406 2407 if (length(extra)) extra$w <- w else extra <- list(w = w) 2408 2409 if (!length(etastart)) { 2410 mu.init <- 2411 if ( .imethod == 1) { 2412 if (length( .iprob )) 2413 matrix( .iprob , n, M, byrow = TRUE) else 2414 if (length(mustart.orig)) 2415 matrix(rep_len(mustart.orig, n * M), n, M) else 2416 mustart # Already n x M 2417 } else { 2418 matrix(runif(n * M), n, M) 2419 } 2420 etastart <- theta2eta(mu.init, .link , earg = .earg ) # n x M 2421 } 2422 mustart <- NULL 2423 }), list( .link = link, .earg = earg, 2424 .type.fitted = type.fitted, 2425 .p.small = p.small, 2426 .no.warning = no.warning, 2427 .iprob = iprob, 2428 .imethod = imethod ))), 2429 2430 linkinv = eval(substitute(function(eta, extra = NULL) { 2431 tau <- extra$ncoly 2432 taup1 <- tau + 1 2433 probs <- eta2theta(eta, .link , earg = .earg ) 2434 prc <- probs[, 1:tau] 2435 prr <- cbind(0, # == pr1.ignored 2436 probs[, taup1:ncol(probs)]) # 1st coln ignored 2437 2438 logQQQ <- rowSums(log1p(-prc)) 2439 QQQ <- exp(logQQQ) 2440 AAA <- exp(log1p(-QQQ)) # 1 - QQQ 2441 2442 type.fitted <- if (length(extra$type.fitted)) extra$type.fitted else { 2443 warning("cannot find 'type.fitted'. ", 2444 "Returning 'likelihood.cond'.") 2445 "likelihood.cond" 2446 } 2447 2448 type.fitted <- match.arg(type.fitted, 2449 c("likelihood.cond", "mean.uncond"))[1] 2450 2451 2452 2453 if ( type.fitted == "likelihood.cond") { 2454 probs.numer <- prr 2455 mat.index <- cbind(1:nrow(prc), extra$cap1) 2456 probs.numer[mat.index] <- prc[mat.index] 2457 probs.numer[extra$cap.hist1 == 0] <- prc[extra$cap.hist1 == 0] 2458 fv <- probs.numer / AAA 2459 } else { 2460 fv <- matrix(prc[, 1] / AAA, nrow(prc), ncol(prc)) 2461 2462 fv[, 2] <- (prc[, 2] + prc[, 1] * (prr[, 2] - prc[, 2])) / AAA 2463 2464 if (tau >= 3) { 2465 QQQcummat <- exp(t( apply(log1p(-prc), 1, cumsum))) 2466 for (jay in 3:tau) { 2467 sum1 <- prc[, 1] 2468 for (kay in 2:(jay-1)) 2469 sum1 <- sum1 + prc[, kay] * QQQcummat[, kay-1] 2470 fv[, jay] <- prc[, jay] * QQQcummat[, jay-1] + 2471 prr[, jay] * sum1 2472 } 2473 fv[, 3:tau] <- fv[, 3:tau] / AAA 2474 } 2475 } 2476 label.cols.y(fv, colnames.y = extra$colnames.y, NOS = NOS) 2477 }, list( .link = link, .earg = earg ))), 2478 last = eval(substitute(expression({ 2479 extra$w <- NULL # Kill it off 2480 2481 2482 misc$link <- rep_len( .link , M) 2483 names(misc$link) <- c(dn2.cap, dn2.recap) 2484 2485 misc$earg <- vector("list", M) 2486 names(misc$earg) <- names(misc$link) 2487 for (ii in 1:M) 2488 misc$earg[[ii]] <- .earg 2489 2490 2491 misc$multiple.responses <- TRUE 2492 misc$iprob <- .iprob 2493 2494 2495 2496 R <- tfit$qr$qr[1:ncol.X.vlm, 1:ncol.X.vlm, drop = FALSE] 2497 R[lower.tri(R)] <- 0 2498 tmp6 <- N.hat.posbernoulli(eta = eta, link = .link , earg = .earg , 2499 R = R, w = w, 2500 X.vlm = X.vlm.save, 2501 Hlist = Hlist, # 20150428; bug fixed here 2502 extra = extra, model.type = "tb") 2503 extra$N.hat <- tmp6$N.hat 2504 extra$SE.N.hat <- tmp6$SE.N.hat 2505 2506 2507 misc$drop.b <- .drop.b 2508 misc$parallel.t <- .parallel.t 2509 misc$parallel.b <- .parallel.b 2510 misc$apply.parint.b <- .apply.parint.b 2511 misc$apply.parint.t <- .apply.parint.t 2512 misc$ridge.constant <- .ridge.constant 2513 misc$ridge.power <- .ridge.power 2514 2515 }), list( .link = link, .earg = earg, 2516 .apply.parint.b = apply.parint.b, 2517 .apply.parint.t = apply.parint.t, 2518 .parallel.t = parallel.t, 2519 .parallel.b = parallel.b, 2520 .drop.b = drop.b, 2521 .ridge.constant = ridge.constant, 2522 .ridge.power = ridge.power, 2523 .iprob = iprob ))), 2524 loglikelihood = eval(substitute( 2525 function(mu, y, w, residuals = FALSE, eta, 2526 extra = NULL, 2527 summation = TRUE) { 2528 2529 tau <- extra$ncoly 2530 taup1 <- tau + 1 2531 ycounts <- y 2532 use.orig.w <- if (length(extra$orig.w)) extra$orig.w else 1 2533 2534 probs <- eta2theta(eta, .link , earg = .earg ) 2535 prc <- probs[, 1:tau] 2536 2537 prr <- cbind(0, # pr1.ignored 2538 probs[, taup1:ncol(probs)]) # 1st coln ignored 2539 2540 2541 if (residuals) { 2542 stop("loglikelihood residuals not implemented yet") 2543 } else { 2544 probs.numer <- prr 2545 mat.index <- cbind(1:nrow(prc), extra$cap1) 2546 probs.numer[mat.index] <- prc[mat.index] 2547 probs.numer[extra$cap.hist1 == 0] <- prc[extra$cap.hist1 == 0] 2548 2549 ll.elts <- 2550 c(use.orig.w) * 2551 dposbern(x = ycounts, # size = 1, # Bernoulli trials 2552 prob = probs.numer, prob0 = prc, log = TRUE) 2553 if (summation) { 2554 sum(ll.elts) 2555 } else { 2556 ll.elts 2557 } 2558 } 2559 }, list( .link = link, .earg = earg ))), 2560 vfamily = c("posbernoulli.tb"), 2561 validparams = eval(substitute(function(eta, y, extra = NULL) { 2562 probs <- eta2theta(eta, .link , earg = .earg ) 2563 okay1 <- all(is.finite(probs)) && all(0 < probs & probs < 1) 2564 okay1 2565 }, list( .link = link, .earg = earg ))), 2566 2567 deriv = eval(substitute(expression({ 2568 tau <- extra$ncoly 2569 taup1 <- tau + 1 2570 probs <- eta2theta(eta, .link , earg = .earg ) 2571 2572 prc <- probs[, 1:tau] 2573 prr <- cbind(pr1.ignored = 0, 2574 probs[, taup1:ncol(probs)]) # 1st coln ignored 2575 logQQQ <- rowSums(log1p(-prc)) 2576 QQQ <- exp(logQQQ) 2577 2578 2579 2580 dprobs.deta <- dtheta.deta(probs, .link , earg = .earg ) 2581 dQ.dprc <- -QQQ / (1 - prc) 2582 d2Q.dprc <- array(0, c(n, tau, tau)) 2583 for (jay in 1:(tau-1)) 2584 for (kay in (jay+1):tau) 2585 d2Q.dprc[, jay, kay] <- 2586 d2Q.dprc[, kay, jay] <- QQQ / ((1 - prc[, jay]) * 2587 (1 - prc[, kay])) 2588 2589 dl.dpc <- dl.dpr <- matrix(0, n, tau) # 1st coln of dl.dpr is ignored 2590 for (jay in 1:tau) { 2591 dl.dpc[, jay] <- (1 - extra$cap.hist1[, jay]) * 2592 ( y[, jay] / prc[, jay] - 2593 (1 - y[, jay]) / (1 - prc[, jay])) + 2594 dQ.dprc[, jay] / (1 - QQQ) 2595 } 2596 for (jay in 2:tau) { 2597 dl.dpr[, jay] <- extra$cap.hist1[, jay] * 2598 ( y[, jay] / prr[, jay] - 2599 (1 - y[, jay]) / (1 - prr[, jay])) 2600 } 2601 2602 deriv.ans <- c(w) * cbind(dl.dpc, dl.dpr[, -1]) * dprobs.deta 2603 deriv.ans 2604 }), list( .link = link, 2605 .earg = earg ))), 2606 2607 weight = eval(substitute(expression({ 2608 wz <- matrix(0, n, sum(M:(M - (tau - 1)))) 2609 2610 2611 2612 QQQcummat <- exp(t( apply(log1p(-prc), 1, cumsum))) 2613 wz.pc <- (QQQcummat / prc - QQQ / (1 - QQQ)) / ((1 - QQQ) * 2614 (1 - prc)^2) 2615 wz[, 1:tau] <- wz.pc 2616 2617 2618 wz.pr <- as.matrix((1 - QQQcummat / (1 - prc)) / ( 2619 prr * (1 - prr) * (1 - QQQ))) 2620 wz[, taup1:M] <- wz.pr[, -1] 2621 2622 2623 for (jay in 1:(tau-1)) 2624 for (kay in (jay+1):tau) 2625 wz[, iam(jay, kay, M = M)] <- 2626 -(d2Q.dprc[, jay, kay] + 2627 dQ.dprc[, jay] * 2628 dQ.dprc[, kay] / (1 - QQQ)) / (1 - QQQ) 2629 2630 2631 cindex <- iam(NA, NA, M = M, both = TRUE) 2632 cindex$row.index <- cindex$row.index[1:ncol(wz)] 2633 cindex$col.index <- cindex$col.index[1:ncol(wz)] 2634 2635 wz <- wz * dprobs.deta[, cindex$row.index] * 2636 dprobs.deta[, cindex$col.index] 2637 2638 2639 if (TRUE) { # ------------------------------------ 2640 wz.adjustment <- .ridge.constant * iter^( .ridge.power ) 2641 wz[, 1:tau] <- wz[, 1:tau] * (1 + wz.adjustment) 2642 } else { # ------------------------------------ 2643 wz.mean <- mean(wz[, 1:tau]) 2644 wz.adjustment <- wz.mean * .ridge.constant * iter^( .ridge.power ) 2645 wz[, 1:tau] <- wz[, 1:tau] + wz.adjustment 2646 } # ------------------------------------ 2647 2648 2649 2650 2651 2652 2653 2654 c(w) * wz 2655 }), list( .link = link, .earg = earg, 2656 .ridge.constant = ridge.constant, 2657 .ridge.power = ridge.power )))) 2658} # posbernoulli.tb 2659 2660 2661 2662 2663 2664 2665setClass("posbernoulli.tb", contains = "vglmff") 2666setClass("posbernoulli.t", contains = "posbernoulli.tb") 2667setClass("posbernoulli.b", contains = "posbernoulli.tb") 2668 2669 setClass("posbinomial", contains = "posbernoulli.b") 2670 2671 2672 2673setMethod("summaryvglmS4VGAM", signature(VGAMff = "posbernoulli.tb"), 2674 function(object, VGAMff, ...) { 2675 object@post 2676}) 2677 2678 2679 2680setMethod("showsummaryvglmS4VGAM", signature(VGAMff = "posbernoulli.tb"), 2681 function(object, VGAMff, ...) { 2682 if (length(object@extra$N.hat) == 1 && 2683 is.numeric(object@extra$N.hat)) { 2684 cat("\nEstimate of N: ", round(object@extra$N.hat, digits = 3), "\n") 2685 cat("\nStd. Error of N: ", round(object@extra$SE.N.hat, digits = 3), 2686 "\n") 2687 2688 confint.N <- object@extra$N.hat + 2689 c(Lower = -1, Upper = 1) * qnorm(0.975) * object@extra$SE.N.hat 2690 cat("\nApproximate 95 percent confidence interval for N:\n") 2691 } 2692}) 2693 2694 2695 2696setMethod("showsummaryvglmS4VGAM", signature(VGAMff = "posbernoulli.b"), 2697 function(object, VGAMff, ...) { 2698 callNextMethod(VGAMff = VGAMff, object = object, ...) 2699}) 2700 2701 2702 2703setMethod("showsummaryvglmS4VGAM", signature(VGAMff = "posbernoulli.t"), 2704 function(object, VGAMff, ...) { 2705 callNextMethod(VGAMff = VGAMff, object = object, ...) 2706}) 2707 2708 2709 2710 2711