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 17hdeffminp <- 18 function(object, ...) { 19 derivs5 <- hdeff(object, deriv = 2, se = FALSE) 20 21 coef(object) - derivs5[, "deriv1"] / derivs5[, "deriv2"] 22} # hdeffminp 23 24 25 26 27 28 29 binomialff <- 30 function(link = "logitlink", 31 multiple.responses = FALSE, 32 parallel = FALSE, # apply.parint = FALSE, 33 zero = NULL, 34 bred = FALSE, 35 earg.link = FALSE) { 36 37 38 dispersion = 1 39 onedpar = !multiple.responses 40 41 42 if (!is.logical(bred) || length(bred) > 1) 43 stop("argument 'bred' must be a single logical") 44 45 46 47 apply.parint <- FALSE 48 estimated.dispersion <- dispersion == 0 49 50 51 52 53 54 if (earg.link) { 55 earg <- link 56 } else { 57 link <- as.list(substitute(link)) 58 earg <- link2list(link) 59 } 60 link <- attr(earg, "function.name") 61 62 63 ans <- 64 new("vglmff", 65 blurb = if (multiple.responses) c("Multiple responses binomial model\n\n", 66 "Link: ", namesof("mu[,j]", link, earg = earg), "\n", 67 "Variance: mu[,j]*(1-mu[,j])") else 68 c("Binomial model\n\n", 69 "Link: ", namesof("prob", link, earg = earg), "\n", 70 "Variance: mu * (1 - mu)"), 71 72 charfun = eval(substitute(function(x, eta, extra = NULL, 73 varfun = FALSE) { 74 mu <- eta2theta(eta, link = .link , earg = .earg ) 75 if (!length(size <- extra$size)) 76 size <- 1 77 if (varfun) { 78 mu * (1 - mu) / size 79 } else { 80 (1 + mu * (exp(1i * x) - 1))^size 81 } 82 }, list( .link = link, .earg = earg ))), 83 84 constraints = eval(substitute(expression({ 85 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 86 bool = .parallel , 87 constraints = constraints, 88 apply.int = .apply.parint ) 89 90 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 91 predictors.names = predictors.names, 92 M1 = 1) 93 }), list( .zero = zero, 94 .parallel = parallel, .apply.parint = apply.parint ))), 95 96 infos = eval(substitute(function(...) { 97 list(M1 = 1, 98 Q1 = 1, 99 bred = .bred , 100 charfun = TRUE, 101 expected = TRUE, 102 hadof = TRUE, 103 multiple.responses = .multiple.responses , 104 parameters.names = c("prob"), # new.name 105 zero = .zero ) 106 }, list( .zero = zero, 107 .bred = bred, 108 .multiple.responses = multiple.responses ))), 109 110 initialize = eval(substitute(expression({ 111 assign("CQO.FastAlgorithm", 112 ( .link == "logitlink" || .link == "clogloglink"), 113 envir = VGAMenv) 114 assign("modelno", if ( .link == "logitlink") 1 else 115 if ( .link == "clogloglink") 4 else NULL, 116 envir = VGAMenv) 117 118 119 120 old.name <- "mu" 121 new.name <- "prob" 122 123 124 125 if ( .multiple.responses ) { 126 temp5 <- 127 w.y.check(w = w, y = y, 128 Is.nonnegative.y = TRUE, 129 ncol.w.max = Inf, 130 ncol.y.max = Inf, 131 out.wy = TRUE, 132 colsyperw = 1, 133 maximize = TRUE) 134 w <- temp5$w 135 y <- temp5$y 136 137 138 y.counts <- y 139 y <- y / w 140 141 142 143 144 145 M <- ncol(y) 146 147 if (FALSE) 148 if (!all(y == 0 | y == 1)) 149 stop("response must contain 0s and 1s only") 150 151 152 153 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 154 dn2 <- if (length(dn2)) { 155 paste("E[", dn2, "]", sep = "") 156 } else { 157 param.names(new.name, M) 158 } 159 predictors.names <- 160 namesof(if (M > 1) dn2 else new.name, 161 .link , earg = .earg , short = TRUE) 162 163 if (!length(mustart) && !length(etastart)) 164 mustart <- matrix(colMeans(y.counts), nrow(y), ncol = ncol(y), 165 byrow = TRUE) / 166 matrix(colMeans(w), nrow = nrow(w), ncol = ncol(w), 167 byrow = TRUE) 168 169 170 171 extra$multiple.responses <- TRUE 172 } else { 173 174 if (!all(w == 1)) 175 extra$orig.w <- w 176 177 178 NCOL <- function (x) if (is.array(x) && length(dim(x)) > 1 || 179 is.data.frame(x)) ncol(x) else as.integer(1) 180 if (NCOL(y) == 1) { 181 if (is.factor(y)) 182 y <- (y != levels(y)[1]) 183 nvec <- rep_len(1, n) 184 y[w == 0] <- 0 185 if (!all(y == 0 | y == 1)) 186 stop("response values 'y' must be 0 or 1") 187 if (!length(mustart) && !length(etastart)) 188 mustart <- (0.5 + w * y) / (1 + w) 189 190 191 no.successes <- y 192 if (min(y) < 0) 193 stop("Negative data not allowed!") 194 if (any(abs(no.successes - round(no.successes)) > 1.0e-8)) 195 stop("Number of successes must be integer-valued") 196 } else if (NCOL(y) == 2) { 197 if (min(y) < 0) 198 stop("Negative data not allowed!") 199 if (any(abs(y - round(y)) > 1.0e-8)) 200 stop("Count data must be integer-valued") 201 y <- round(y) 202 nvec <- y[, 1] + y[, 2] 203 y <- ifelse(nvec > 0, y[, 1] / nvec, 0) 204 w <- w * nvec 205 if (!length(mustart) && !length(etastart)) 206 mustart <- (0.5 + nvec * y) / (1 + nvec) 207 } else { 208 stop("for the binomialff family, response 'y' must be a ", 209 "vector of 0 and 1's\n", 210 "or a factor (first level = fail, other levels = success)", 211 ",\n", 212 "or a 2-column matrix where col 1 is the no. of ", 213 "successes and col 2 is the no. of failures") 214 } 215 predictors.names <- 216 namesof(new.name, .link , earg = .earg , short = TRUE) 217 } # Not multiple.responses. 218 219 220 if ( .bred ) { 221 if ( !control$save.weights ) { 222 save.weights <- control$save.weights <- TRUE 223 } 224 } 225 226 227 228 }), list( .link = link, .multiple.responses = multiple.responses, 229 .earg = earg, .bred = bred ))), 230 231 linkinv = eval(substitute(function(eta, extra = NULL) { 232 mu <- eta2theta(eta, link = .link , earg = .earg ) 233 234 235 colnames(mu) <- NULL 236 237 mu 238 }, list( .link = link, .earg = earg ))), 239 240 last = eval(substitute(expression({ 241 if (exists("CQO.FastAlgorithm", envir = VGAMenv)) 242 rm("CQO.FastAlgorithm", envir = VGAMenv) 243 if (exists("modelno", envir = VGAMenv)) 244 rm("modelno", envir = VGAMenv) 245 246 dpar <- .dispersion 247 if (!dpar) { 248 temp87 <- (y-mu)^2 * wz / ( 249 dtheta.deta(mu, link = .link , 250 earg = .earg )^2) # w cancel 251 if (.multiple.responses && ! .onedpar ) { 252 dpar <- rep_len(NA_real_, M) 253 temp87 <- cbind(temp87) 254 nrow.mu <- if (is.matrix(mu)) nrow(mu) else length(mu) 255 for (ii in 1:M) 256 dpar[ii] <- sum(temp87[, ii]) / (nrow.mu - ncol(x)) 257 if (is.matrix(y) && length(dimnames(y)[[2]]) == length(dpar)) 258 names(dpar) <- dimnames(y)[[2]] 259 } else { 260 dpar <- sum(temp87) / (length(mu) - ncol(x)) 261 } 262 } 263 264 265 266 if (! ( .multiple.responses )) 267 extra$size <- nvec # For @charfun, and therefore "stdres" 268 269 misc$dispersion <- dpar 270 misc$default.dispersion <- 1 271 misc$estimated.dispersion <- .estimated.dispersion 272 273 misc$link <- rep_len( .link , M) 274 names(misc$link) <- if (M > 1) dn2 else new.name # Was old.name=="mu" 275 276 misc$earg <- vector("list", M) 277 names(misc$earg) <- names(misc$link) 278 for (ii in 1:M) 279 misc$earg[[ii]] <- .earg 280 281 }), list( .dispersion = dispersion, 282 .estimated.dispersion = estimated.dispersion, 283 .onedpar = onedpar, .multiple.responses = multiple.responses, 284 .bred = bred, 285 .link = link, .earg = earg))), 286 287 linkfun = eval(substitute(function(mu, extra = NULL) { 288 theta2eta(mu, .link , earg = .earg ) 289 }, list( .link = link, .earg = earg))), 290 291 loglikelihood = eval(substitute( 292 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 293 summation = TRUE) { 294 if (residuals) { 295 c(w) * (y / mu - (1-y) / (1-mu)) 296 } else { 297 298 ycounts <- if (is.numeric(extra$orig.w)) y * w / extra$orig.w else 299 y * w # Convert proportions to counts 300 nvec <- if (is.numeric(extra$orig.w)) round(w / extra$orig.w) else 301 round(w) 302 303 304 smallno <- 1.0e6 * .Machine$double.eps 305 smallno <- sqrt(.Machine$double.eps) 306 if (max(abs(ycounts - round(ycounts))) > smallno) 307 warning("converting 'ycounts' to integer in @loglikelihood") 308 ycounts <- round(ycounts) 309 310 ll.elts <- if ( .multiple.responses ) { 311 c(w) * ( ycounts * log( mu) + 312 (1 - ycounts) * log1p(-mu)) 313 } else { 314 (if (is.numeric(extra$orig.w)) extra$orig.w else 1) * 315 dbinom(x = ycounts, size = nvec, prob = mu, log = TRUE) 316 } 317 if (summation) { 318 sum(ll.elts) 319 } else { 320 ll.elts 321 } 322 } 323 }, list( .multiple.responses = multiple.responses ))), 324 325 validparams = eval(substitute(function(eta, y, extra = NULL) { 326 mymu <- eta2theta(eta, .link , earg = .earg ) 327 okay1 <- all(is.finite(mymu)) && all(0 < mymu & mymu < 1) 328 okay1 329 }, list( .link = link, .earg = earg, .bred = bred))), 330 vfamily = c("binomialff", "VGAMglm"), # "VGAMcategorical", 331 332 333 334 335 hadof = eval(substitute( 336 function(eta, extra = list(), deriv = 1, 337 linpred.index = 1, 338 w = 1, dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2), 339 ...) { 340 fvs <- eta2theta(eta, link = .link , earg = .earg ) 341 if ( .bred ) { 342 fvs <- fvs + NA_real_ # Correct dimension for below too 343 } 344 345 ans <- c(w) * 346 switch(as.character(deriv), 347 "0" = 1 / (fvs * (1 - fvs)), 348 "1" = -(1 - 2*fvs) / (fvs * (1 - fvs))^2, 349 "2" = 2 * (1 - 3*fvs*(1-fvs)) / (fvs * (1 - fvs))^3, 350 stop("argument 'deriv' must be 0 or 1 or 2")) 351 if (deriv == 0) ans else retain.col(ans, linpred.index) # Coz M1 = 1 352 }, list( .link = link, .earg = earg, .bred = bred) )), 353 354 355 356 357 358 simslot = function (object, nsim) { 359 360 ftd <- fitted(object) 361 362 363 if (ncol(ftd) > 1) 364 stop("simulate() does not work with more than one response") 365 366 367 368 n <- length(ftd) 369 ntot <- n * nsim 370 pwts <- if (length(pwts <- object@prior.weights) > 0) 371 pwts else weights(object, type = "prior") 372 373 374 375 if (any(pwts %% 1 != 0)) 376 stop("cannot simulate from non-integer prior.weights") 377 if (length(m <- object@model) > 0) { 378 y <- model.response(m) 379 if (is.factor(y)) { 380 yy <- factor(1 + rbinom(ntot, size = 1, prob = ftd), 381 labels = levels(y)) 382 split(yy, rep(seq_len(nsim), each = n)) 383 } else if (is.matrix(y) && ncol(y) == 2) { 384 yy <- vector("list", nsim) 385 for (i in seq_len(nsim)) { 386 Y <- rbinom(n, size = pwts, prob = ftd) 387 YY <- cbind(Y, pwts - Y) 388 colnames(YY) <- colnames(y) 389 yy[[i]] <- YY 390 } 391 yy 392 } else { 393 rbinom(ntot, size = pwts, prob = ftd)/pwts 394 } 395 } else { 396 397 rbinom(ntot, size = c(pwts), prob = c(ftd))/c(pwts) 398 } 399 }, 400 401 402 403 404 deriv = eval(substitute(expression({ 405 406 407 408 409 410 411 yBRED <- if ( .bred ) { 412 Hvector <- hatvaluesbasic(X.vlm = X.vlm.save, 413 diagWm = c(t(w * mu))) # Handles M>1 414 415 varY <- mu * (1 - mu) / w # A matrix if M>1. Seems the most correct. 416 d1.ADJ <- dtheta.deta(mu, .link , earg = .earg ) 417 418 temp.earg <- .earg 419 temp.earg$inverse <- FALSE 420 temp.earg$inverse <- TRUE 421 d2.ADJ <- d2theta.deta2(mu, .link , earg = temp.earg ) 422 423 424 425 yBRED <- y + matrix(Hvector, n, M, byrow = TRUE) * 426 varY * d2.ADJ / (2 * d1.ADJ^2) 427 yBRED 428 } else { 429 y 430 } 431 432 433 answer <- if ( .link == "logitlink") { 434 c(w) * (yBRED - mu) 435 } else if ( .link == "clogloglink") { 436 mu.use <- mu 437 smallno <- 100 * .Machine$double.eps 438 mu.use[mu.use < smallno] <- smallno 439 mu.use[mu.use > 1.0 - smallno] <- 1.0 - smallno 440 -c(w) * (yBRED - mu) * log1p(-mu.use) / mu.use 441 } else { 442 c(w) * dtheta.deta(mu, link = .link , earg = .earg ) * 443 (yBRED / mu - 1.0) / (1.0 - mu) 444 } 445 446 answer 447 }), list( .link = link, .earg = earg, .bred = bred))), 448 449 weight = eval(substitute(expression({ 450 tmp100 <- mu * (1.0 - mu) 451 452 ned2ldprob2 <- if ( .link == "logitlink") { 453 cbind(c(w) * tmp100) 454 } else if ( .link == "clogloglink") { 455 cbind(c(w) * (1.0 - mu.use) * (log1p(-mu.use))^2 / mu.use) 456 } else { 457 cbind(c(w) * dtheta.deta(mu, .link , earg = .earg )^2 / tmp100) 458 } 459 for (ii in 1:M) { 460 index500 <- !is.finite(ned2ldprob2[, ii]) | 461 (abs(ned2ldprob2[, ii]) < .Machine$double.eps) 462 if (any(index500)) { # Diagonal 0s are bad 463 ned2ldprob2[index500, ii] <- .Machine$double.eps 464 } 465 } 466 ned2ldprob2 467 }), list( .link = link, .earg = earg)))) 468 469 470 471 472 473 474 ans@deviance <- 475 if (multiple.responses) 476 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 477 summation = TRUE) { 478 479 M <- NOS <- NCOL(mu) 480 mu.use <- cbind(mu, 1 - mu)[, interleave.VGAM(2 * M, M1 = 2), 481 drop = FALSE] 482 yy.use <- cbind( y, 1 - y)[, interleave.VGAM(2 * M, M1 = 2), 483 drop = FALSE] 484 ww.use <- kronecker(matrix(w, NROW(y), M), cbind(1, 1)) 485 Deviance.categorical.data.vgam(mu = mu.use, 486 y = yy.use, 487 w = ww.use, residuals = residuals, 488 eta = eta, extra = extra, 489 summation = summation) 490 } else 491 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 492 summation = TRUE) { 493 Deviance.categorical.data.vgam(mu = cbind(mu, 1-mu), 494 y = cbind(y , 1-y), 495 w = w, residuals = residuals, 496 eta = eta, extra = extra, 497 summation = summation) 498 } 499 ans 500} # binomialff() 501 502 503 504 505 506 gammaff <- function(link = "negreciprocal", dispersion = 0) { 507 estimated.dispersion <- dispersion == 0 508 509 510 link <- as.list(substitute(link)) 511 earg <- link2list(link) 512 link <- attr(earg, "function.name") 513 514 515 new("vglmff", 516 blurb = c("Gamma distribution\n\n", 517 "Link: ", namesof("mu", link, earg = earg), "\n", 518 "Variance: mu^2 / k"), 519 deviance = 520 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 521 summation = TRUE) { 522 devi <- -2 * c(w) * (log(ifelse(y == 0, 1, y/mu)) - (y - mu)/mu) 523 if (residuals) { 524 sign(y - mu) * sqrt(abs(devi) * w) 525 } else { 526 dev.elts <- c(w) * devi 527 if (summation) { 528 sum(dev.elts) 529 } else { 530 dev.elts 531 } 532 } 533 }, 534 infos = eval(substitute(function(...) { 535 list(M1 = 1, 536 Q1 = 1, 537 parameters.names = c("mu"), 538 dispersion = .dispersion ) 539 }, list( .dispersion = dispersion ))), 540 initialize = eval(substitute(expression({ 541 542 temp5 <- 543 w.y.check(w = w, y = y, 544 Is.nonnegative.y = TRUE, 545 out.wy = TRUE, 546 colsyperw = 1, 547 maximize = TRUE) 548 w <- temp5$w 549 y <- temp5$y 550 551 552 mustart <- y + 0.167 * (y == 0) 553 554 M <- if (is.matrix(y)) ncol(y) else 1 555 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 556 dn2 <- if (length(dn2)) { 557 paste("E[", dn2, "]", sep = "") 558 } else { 559 param.names("mu", M, skip1 = TRUE) 560 } 561 562 predictors.names <- 563 namesof(if (M > 1) dn2 else "mu", .link , 564 earg = .earg , short = TRUE) 565 566 if (!length(etastart)) 567 etastart <- theta2eta(mustart, link = .link , earg = .earg ) 568 }), list( .link = link, .earg = earg))), 569 linkinv = eval(substitute(function(eta, extra = NULL) { 570 eta2theta(eta, link = .link , earg = .earg ) 571 }, list( .link = link, .earg = earg))), 572 last = eval(substitute(expression({ 573 dpar <- .dispersion 574 if (!dpar) { 575 if (M == 1) { 576 temp <- c(w) * dmu.deta^2 577 dpar <- sum(c(w) * (y-mu)^2 * wz / temp) / (length(mu) - ncol(x)) 578 } else { 579 dpar <- rep_len(0, M) 580 for (spp in 1:M) { 581 temp <- c(w) * dmu.deta[, spp]^2 582 dpar[spp] <- sum(c(w) * (y[,spp]-mu[, spp])^2 * 583 wz[, spp]/temp) / ( 584 length(mu[,spp]) - ncol(x)) 585 } 586 } 587 } 588 misc$dispersion <- dpar 589 misc$default.dispersion <- 0 590 misc$estimated.dispersion <- .estimated.dispersion 591 592 misc$link <- rep_len( .link , M) 593 names(misc$link) <- param.names("mu", M, skip1 = TRUE) 594 595 misc$earg <- vector("list", M) 596 names(misc$earg) <- names(misc$link) 597 for (ii in 1:M) 598 misc$earg[[ii]] <- .earg 599 600 misc$expected <- TRUE 601 misc$multipleResponses <- TRUE 602 }), list( .dispersion = dispersion, .earg = earg, 603 .estimated.dispersion = estimated.dispersion, 604 .link = link ))), 605 linkfun = eval(substitute(function(mu, extra = NULL) { 606 theta2eta(mu, link = .link , earg = .earg ) 607 }, list( .link = link, .earg = earg))), 608 vfamily = "gammaff", 609 validparams = eval(substitute(function(eta, y, extra = NULL) { 610 mymu <- theta2eta(mu, .link , earg = .earg ) 611 okay1 <- all(is.finite(mymu)) && all(0 < mymu) 612 okay1 613 }, list( .link = link, .earg = earg ))), 614 deriv = eval(substitute(expression({ 615 M1 <- 1 616 ncoly <- NCOL(y) 617 618 dl.dmu <- (y-mu) / mu^2 619 dmu.deta <- dtheta.deta(theta = mu, link = .link , earg = .earg ) 620 c(w) * dl.dmu * dmu.deta 621 }), list( .link = link, .earg = earg))), 622 weight = eval(substitute(expression({ 623 d2l.dmu2 <- 1 / mu^2 624 wz <- dmu.deta^2 * d2l.dmu2 625 w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = ncoly) 626 }), list( .link = link, .earg = earg)))) 627} # gammaff() 628 629 630 631 632 633 634 inverse.gaussianff <- function(link = "natural.ig", 635 dispersion = 0) { 636 estimated.dispersion <- dispersion == 0 637 warning("@deviance() not finished") 638 warning("needs checking, but I'm sure it works") 639 640 link <- as.list(substitute(link)) 641 earg <- link2list(link) 642 link <- attr(earg, "function.name") 643 644 645 new("vglmff", 646 blurb = c("Inverse Gaussian distribution\n\n", 647 "Link: ", namesof("mu", link, earg = earg), "\n", 648 "Variance: mu^3 / k"), 649 650 deviance = 651 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 652 summation = TRUE) { 653 pow <- 3 # Use Quasi()$deviance with pow==3 654 devy <- y^(2-pow) / (1-pow) - y^(2-pow) / (2-pow) 655 devmu <- y * mu^(1-pow) / (1-pow) - mu^(2-pow) / (2-pow) 656 devi <- 2 * (devy - devmu) 657 if (residuals) { 658 sign(y - mu) * sqrt(abs(devi) * w) 659 } else { 660 dev.elts <- c(w) * devi 661 if (summation) { 662 sum(dev.elts) 663 } else { 664 dev.elts 665 } 666 } 667 }, 668 669 infos = eval(substitute(function(...) { 670 list(M1 = 1, 671 Q1 = 1, 672 parameters.names = c("mu"), 673 quasi.type = TRUE, 674 dispersion = .dispersion ) 675 }, list( .earg = earg , .dispersion = dispersion ))), 676 initialize = eval(substitute(expression({ 677 temp5 <- 678 w.y.check(w = w, y = y, 679 Is.positive.y = TRUE, 680 out.wy = TRUE, 681 colsyperw = 1, 682 maximize = TRUE) 683 w <- temp5$w 684 y <- temp5$y 685 686 687 mu <- y + 0.167 * (y == 0) 688 689 690 691 M <- if (is.matrix(y)) ncol(y) else 1 692 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 693 dn2 <- if (length(dn2)) { 694 paste("E[", dn2, "]", sep = "") 695 } else { 696 param.names("mu", M, skip1 = TRUE) 697 } 698 699 predictors.names <- 700 namesof(if (M > 1) dn2 else "mu", .link , .earg , short = TRUE) 701 702 703 if (!length(etastart)) 704 etastart <- theta2eta(mu, link = .link , .earg ) 705 }), list( .link = link, .earg = earg))), 706 linkinv = eval(substitute(function(eta, extra = NULL) { 707 eta2theta(eta, link = .link , earg = .earg ) 708 }, list( .link = link, .earg = earg ))), 709 last = eval(substitute(expression({ 710 dpar <- .dispersion 711 if (!dpar) { 712 temp <- c(w) * dmu.deta^2 713 dpar <- sum( c(w) * (y-mu)^2 * wz / temp ) / (length(mu) - ncol(x)) 714 } 715 misc$dispersion <- dpar 716 misc$default.dispersion <- 0 717 misc$estimated.dispersion <- .estimated.dispersion 718 719 misc$link <- rep_len( .link , M) 720 names(misc$link) <- param.names("mu", M, skip1 = TRUE) 721 722 misc$earg <- vector("list", M) 723 names(misc$earg) <- names(misc$link) 724 for (ii in 1:M) 725 misc$earg[[ii]] <- .earg 726 727 misc$expected <- TRUE 728 misc$multipleResponses <- TRUE 729 }), list( .dispersion = dispersion, 730 .estimated.dispersion = estimated.dispersion, 731 .link = link, .earg = earg ))), 732 linkfun = eval(substitute(function(mu, extra = NULL) { 733 theta2eta(mu, link = .link , earg = .earg ) 734 }, list( .link = link, .earg = earg ))), 735 vfamily = "inverse.gaussianff", 736 validparams = eval(substitute(function(eta, y, extra = NULL) { 737 mymu <- theta2eta(mu, .link , earg = .earg ) 738 okay1 <- all(is.finite(mymu)) && all(0 < mymu) 739 okay1 740 }, list( .link = link, .earg = earg ))), 741 deriv = eval(substitute(expression({ 742 M1 <- 1 743 ncoly <- NCOL(y) 744 745 dl.dmu <- (y - mu) / mu^3 746 dmu.deta <- dtheta.deta(mu, link = .link , earg = .earg ) 747 c(w) * dl.dmu * dmu.deta 748 }), list( .link = link, .earg = earg ))), 749 weight = eval(substitute(expression({ 750 d2l.dmu2 <- 1 / mu^3 751 wz <- dmu.deta^2 * d2l.dmu2 752 w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = ncoly) 753 }), list( .link = link, .earg = earg )))) 754} 755 756 757 758 759dinv.gaussian <- function(x, mu, lambda, log = FALSE) { 760 if (!is.logical(log.arg <- log) || length(log) != 1) 761 stop("bad input for argument 'log'") 762 rm(log) 763 764 L <- max(length(x), length(mu), length(lambda)) 765 if (length(x) != L) x <- rep_len(x, L) 766 if (length(mu) != L) mu <- rep_len(mu, L) 767 if (length(lambda) != L) lambda <- rep_len(lambda, L) 768 logdensity <- rep_len(log(0), L) 769 770 xok <- (x > 0) 771 logdensity[xok] <- 0.5 * log(lambda[xok] / (2 * pi * x[xok]^3)) - 772 lambda[xok] * (x[xok]-mu[xok])^2 / (2*mu[xok]^2 * x[xok]) 773 logdensity[mu <= 0] <- NaN 774 logdensity[lambda <= 0] <- NaN 775 if (log.arg) logdensity else exp(logdensity) 776} 777 778 779 780pinv.gaussian <- function(q, mu, lambda) { 781 L <- max(length(q), length(mu), length(lambda)) 782 if (length(q) != L) q <- rep_len(q, L) 783 if (length(mu) != L) mu <- rep_len(mu, L) 784 if (length(lambda) != L) lambda <- rep_len(lambda, L) 785 ans <- q 786 787 ans[q <= 0] <- 0 788 bb <- q > 0 789 ans[bb] <- pnorm( sqrt(lambda[bb]/q[bb]) * (q[bb]/mu[bb] - 1)) + 790 exp(2*lambda[bb]/mu[bb]) * 791 pnorm(-sqrt(lambda[bb]/q[bb]) * (q[bb]/mu[bb] + 1)) 792 ans[mu <= 0] <- NaN 793 ans[lambda <= 0] <- NaN 794 ans 795} 796 797 798 799rinv.gaussian <- function(n, mu, lambda) { 800 use.n <- if ((length.n <- length(n)) > 1) length.n else 801 if (!is.Numeric(n, integer.valued = TRUE, 802 length.arg = 1, positive = TRUE)) 803 stop("bad input for argument 'n'") else n 804 805 mu <- rep_len(mu, use.n) 806 lambda <- rep_len(lambda, use.n) 807 808 u <- runif(use.n) 809 Z <- rnorm(use.n)^2 # rchisq(use.n, df = 1) 810 phi <- lambda / mu 811 y1 <- 1 - 0.5 * (sqrt(Z^2 + 4*phi*Z) - Z) / phi 812 ans <- mu * ifelse((1+y1)*u > 1, 1/y1, y1) 813 ans[mu <= 0] <- NaN 814 ans[lambda <= 0] <- NaN 815 ans 816} 817 818 819 820 821 822 823 824 inv.gaussianff <- function(lmu = "loglink", llambda = "loglink", 825 imethod = 1, ilambda = NULL, 826 parallel = FALSE, 827 ishrinkage = 0.99, 828 zero = NULL) { 829 830 831 apply.parint <- FALSE 832 833 834 lmu <- as.list(substitute(lmu)) 835 emu <- link2list(lmu) 836 lmu <- attr(emu, "function.name") 837 838 llambda <- as.list(substitute(llambda)) 839 elambda <- link2list(llambda) 840 llambda <- attr(elambda, "function.name") 841 842 843 if (!is.Numeric(imethod, length.arg = 1, 844 integer.valued = TRUE, positive = TRUE) || 845 imethod > 3) 846 stop("argument 'imethod' must be 1 or 2 or 3") 847 if (!is.Numeric(ishrinkage, length.arg = 1) || 848 ishrinkage < 0 || 849 ishrinkage > 1) 850 stop("bad input for argument 'ishrinkage'") 851 852 853 if (is.logical(parallel) && parallel && length(zero)) 854 stop("set 'zero = NULL' if 'parallel = TRUE'") 855 856 857 858 859 new("vglmff", 860 blurb = c("Inverse Gaussian distribution\n\n", 861 "f(y) = sqrt(lambda/(2*pi*y^3)) * ", 862 "exp(-lambda * (y - mu)^2 / (2 * mu^2 * y)); ", 863 "y, mu & lambda > 0", 864 "Link: ", namesof("mu", lmu, earg = emu), ", ", 865 namesof("lambda", llambda, earg = elambda),"\n", 866 "Mean: ", "mu\n", 867 "Variance: mu^3 / lambda"), 868 constraints = eval(substitute(expression({ 869 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 870 bool = .parallel , 871 constraints = constraints, 872 apply.int = .apply.parint ) 873 874 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 875 predictors.names = predictors.names, 876 M1 = 2) 877 }), list( .zero = zero, 878 .parallel = parallel, .apply.parint = apply.parint ))), 879 infos = eval(substitute(function(...) { 880 list(M1 = 2, 881 Q1 = 1, 882 imethod = .imethod , 883 parameters.names = c("mu", "lambda"), 884 expected = TRUE, 885 multipleResponses = FALSE, 886 zero = .zero ) 887 }, list( .imethod = imethod, .zero = zero ))), 888 889 initialize = eval(substitute(expression({ 890 temp5 <- 891 w.y.check(w = w, y = y, 892 Is.positive.y = TRUE, 893 ncol.w.max = Inf, 894 ncol.y.max = Inf, 895 out.wy = TRUE, 896 colsyperw = 1, 897 maximize = TRUE) 898 w <- temp5$w 899 y <- temp5$y 900 901 902 ncoly <- ncol(y) 903 M1 <- 2 904 extra$ncoly <- ncoly 905 extra$M1 <- M1 906 M <- M1 * ncoly 907 908 909 910 mynames1 <- param.names("mu", ncoly, skip1 = TRUE) 911 mynames2 <- param.names("lambda", ncoly, skip1 = TRUE) 912 predictors.names <- 913 c(namesof(mynames1, .lmu , earg = .emu , short = TRUE), 914 namesof(mynames2, .llambda , earg = .elambda , short = TRUE))[ 915 interleave.VGAM(M, M1 = M1)] 916 917 918 919 920 if (!length(etastart)) { 921 init.mu <- 922 if ( .imethod == 2) { 923 mediany <- apply(y, 2, median) 924 matrix(1.1 * mediany + 1/8, n, ncoly, byrow = TRUE) 925 } else if ( .imethod == 3) { 926 use.this <- colSums(y * w) / colSums(w) # weighted.mean(y, w) 927 (1 - .ishrinkage ) * y + .ishrinkage * use.this 928 } else { 929 matrix(colSums(y * w) / colSums(w) + 1/8, 930 n, ncoly, byrow = TRUE) 931 } 932 933 variancey <- apply(y, 2, var) 934 init.la <- matrix(if (length( .ilambda )) .ilambda else 935 (init.mu^3) / (0.10 + variancey), 936 n, ncoly, byrow = TRUE) 937 938 etastart <- cbind( 939 theta2eta(init.mu, link = .lmu , earg = .emu ), 940 theta2eta(init.la, link = .llambda , earg = .elambda ))[, 941 interleave.VGAM(M, M1 = M1)] 942 } 943 }), list( .lmu = lmu, .llambda = llambda, 944 .emu = emu, .elambda = elambda, 945 .ishrinkage = ishrinkage, 946 .imethod = imethod, .ilambda = ilambda ))), 947 948 linkinv = eval(substitute(function(eta, extra = NULL) { 949 eta2theta(eta[, c(TRUE, FALSE)], link = .lmu , earg = .emu ) 950 }, list( .lmu = lmu, .emu = emu, .elambda = elambda ))), 951 952 last = eval(substitute(expression({ 953 M1 <- extra$M1 954 misc$link <- c(rep_len( .lmu , ncoly), 955 rep_len( .llambda , ncoly))[interleave.VGAM(M, M1 = M1)] 956 temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)] 957 names(misc$link) <- temp.names 958 959 misc$earg <- vector("list", M) 960 names(misc$earg) <- temp.names 961 for (ii in 1:ncoly) { 962 misc$earg[[M1*ii-1]] <- .emu 963 misc$earg[[M1*ii ]] <- .elambda 964 } 965 966 misc$ishrinkage <- .ishrinkage 967 misc$parallel <- .parallel 968 misc$apply.parint <- .apply.parint 969 }), list( .lmu = lmu, .llambda = llambda, 970 .emu = emu, .elambda = elambda, 971 .parallel = parallel, .apply.parint = apply.parint, 972 .ishrinkage = ishrinkage, 973 .imethod = imethod ))), 974 975 loglikelihood = eval(substitute( 976 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 977 summation = TRUE) { 978 mymu <- eta2theta(eta[, c(TRUE, FALSE)], 979 link = .lmu , earg = .emu ) 980 lambda <- eta2theta(eta[, c(FALSE, TRUE)], 981 link = .llambda , earg = .elambda ) 982 if (residuals) { 983 stop("loglikelihood residuals not implemented yet") 984 } else { 985 ll.elts <- c(w) * dinv.gaussian(x = y, mu = mymu, 986 lambda = lambda, log = TRUE) 987 if (summation) { 988 sum(ll.elts) 989 } else { 990 ll.elts 991 } 992 } 993 }, list( .lmu = lmu, .llambda = llambda, 994 .emu = emu, .elambda = elambda ))), 995 996 vfamily = "inv.gaussianff", 997 validparams = eval(substitute(function(eta, y, extra = NULL) { 998 mymu <- eta2theta(eta[, c(TRUE, FALSE)], 999 link = .lmu , earg = .emu ) 1000 lambda <- eta2theta(eta[, c(FALSE, TRUE)], 1001 link = .llambda , earg = .elambda ) 1002 okay1 <- all(is.finite(mymu )) && all(0 < mymu ) && 1003 all(is.finite(lambda)) && all(0 < lambda) 1004 okay1 1005 }, list( .lmu = lmu, .llambda = llambda, 1006 .emu = emu, .elambda = elambda ))), 1007 1008 deriv = eval(substitute(expression({ 1009 M1 <- 2 1010 mymu <- eta2theta(eta[, c(TRUE, FALSE)], 1011 link = .lmu , earg = .emu ) 1012 lambda <- eta2theta(eta[, c(FALSE, TRUE)], 1013 link = .llambda , earg = .elambda ) 1014 1015 dmu.deta <- dtheta.deta(theta = mymu , link = .lmu , earg = .emu ) 1016 dlambda.deta <- dtheta.deta(theta = lambda, link = .llambda , 1017 earg = .elambda ) 1018 1019 dl.dmu <- lambda * (y - mymu) / mymu^3 1020 dl.dlambda <- 0.5 / lambda - (y - mymu)^2 / (2 * mymu^2 * y) 1021 myderiv <- c(w) * cbind(dl.dmu * dmu.deta, 1022 dl.dlambda * dlambda.deta) 1023 myderiv[, interleave.VGAM(M, M1 = M1)] 1024 }), list( .lmu = lmu, .llambda = llambda, 1025 .emu = emu, .elambda = elambda ))), 1026 1027 weight = eval(substitute(expression({ 1028 1029 ned2l.dmu2 <- lambda / mymu^3 1030 ned2l.dlambda2 <- 0.5 / (lambda^2) 1031 1032 wz <- cbind(dmu.deta^2 * ned2l.dmu2, 1033 dlambda.deta^2 * ned2l.dlambda2)[, 1034 interleave.VGAM(M, M1 = M1)] 1035 1036 w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1) 1037 }), list( .lmu = lmu, .llambda = llambda, 1038 .emu = emu, .elambda = elambda )))) 1039} # inv.gaussianff() 1040 1041 1042 1043 1044 1045 poissonff <- function(link = "loglink", 1046 imu = NULL, imethod = 1, 1047 parallel = FALSE, zero = NULL, 1048 bred = FALSE, 1049 earg.link = FALSE, 1050 type.fitted = c("mean", "quantiles"), 1051 percentiles = c(25, 50, 75)) { 1052 1053 1054 1055 dispersion = 1 1056 onedpar = FALSE 1057 1058 1059 type.fitted <- match.arg(type.fitted, 1060 c("mean", "quantiles"))[1] 1061 1062 if (!is.logical(bred) || length(bred) > 1) 1063 stop("argument 'bred' must be a single logical") 1064 1065 estimated.dispersion <- (dispersion == 0) 1066 1067 1068 if (earg.link) { 1069 earg <- link 1070 } else { 1071 link <- as.list(substitute(link)) 1072 earg <- link2list(link) 1073 } 1074 link <- attr(earg, "function.name") 1075 1076 1077 1078 1079 if (!is.Numeric(imethod, length.arg = 1, 1080 integer.valued = TRUE, positive = TRUE) || 1081 imethod > 3) 1082 stop("argument 'imethod' must be 1 or 2 or 3") 1083 if (length(imu) && 1084 !is.Numeric(imu, positive = TRUE)) 1085 stop("bad input for argument 'imu'") 1086 1087 1088 new("vglmff", 1089 blurb = c("Poisson distribution\n\n", 1090 "Link: ", namesof("lambda", link, earg = earg), "\n", 1091 "Variance: lambda"), 1092 1093 charfun = eval(substitute(function(x, eta, extra = NULL, 1094 varfun = FALSE) { 1095 lambda <- eta2theta(eta, link = .link , earg = .earg ) 1096 if (varfun) { 1097 lambda 1098 } else { 1099 exp(lambda * (exp(1i * x) - 1)) 1100 } 1101 }, list( .link = link, .earg = earg ))), 1102 1103 constraints = eval(substitute(expression({ 1104 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1105 bool = .parallel , 1106 constraints = constraints) 1107 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1108 predictors.names = predictors.names, 1109 M1 = 1) 1110 }), list( .parallel = parallel, .zero = zero ))), 1111 1112 infos = eval(substitute(function(...) { 1113 list(M1 = 1, 1114 Q1 = 1, 1115 charfun = TRUE, 1116 expected = TRUE, 1117 hadof = TRUE, 1118 multipleResponses = TRUE, 1119 parameters.names = c("lambda"), 1120 type.fitted = .type.fitted , 1121 percentiles = .percentiles , 1122 bred = .bred , 1123 zero = .zero ) 1124 }, list( .zero = zero, 1125 .type.fitted = type.fitted, 1126 .percentiles = percentiles, 1127 .bred = bred ))), 1128 1129 1130 deviance = eval(substitute( 1131 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 1132 summation = TRUE) { 1133 mupo <- eta2theta(eta, link = .link , earg = .earg ) 1134 nz <- (y > 0) 1135 devi <- -(y - mupo) 1136 devi[nz] <- devi[nz] + y[nz] * log(y[nz]/mupo[nz]) 1137 if (residuals) { 1138 sign(y - mupo) * sqrt(2 * abs(devi) * c(w)) 1139 } else { 1140 dev.elts <- 2 * c(w) * devi 1141 if (summation) { 1142 sum(dev.elts) 1143 } else { 1144 dev.elts 1145 } 1146 } 1147 }, list( .link = link, .earg = earg ))), 1148 1149 initialize = eval(substitute(expression({ 1150 1151 temp5 <- 1152 w.y.check(w = w, y = y, 1153 Is.nonnegative.y = TRUE, 1154 ncol.w.max = Inf, 1155 ncol.y.max = Inf, 1156 out.wy = TRUE, 1157 colsyperw = 1, 1158 maximize = TRUE) 1159 w <- temp5$w 1160 y <- temp5$y 1161 1162 1163 M <- ncoly <- ncol(y) 1164 1165 assign("CQO.FastAlgorithm", ( .link == "loglink"), envir = VGAMenv) 1166 1167 1168 1169 old.name <- "mu" 1170 new.name <- "lambda" 1171 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 1172 dn2 <- if (length(dn2)) { 1173 paste("E[", dn2, "]", sep = "") 1174 } else { 1175 param.names(new.name, M) 1176 } 1177 predictors.names <- 1178 namesof(if (M > 1) dn2 else new.name, # was "mu" == old.name 1179 .link , 1180 earg = .earg , short = TRUE) 1181 1182 1183 if ( .bred ) { 1184 if ( !control$save.weights ) { 1185 save.weights <- control$save.weights <- TRUE 1186 } 1187 } 1188 1189 extra$type.fitted <- .type.fitted 1190 extra$colnames.y <- colnames(y) 1191 extra$percentiles <- .percentiles 1192 1193 1194 if (!length(etastart)) { 1195 mu.init <- pmax(y, 1/8) 1196 for (iii in 1:ncol(y)) { 1197 if ( .imethod == 2) { 1198 mu.init[, iii] <- weighted.mean(y[, iii], w[, iii]) + 1/8 1199 } else if ( .imethod == 3) { 1200 mu.init[, iii] <- median(y[, iii]) + 1/8 1201 } 1202 } 1203 if (length( .imu )) 1204 mu.init <- matrix( .imu , n, ncoly, byrow = TRUE) 1205 etastart <- theta2eta(mu.init, link = .link , earg = .earg ) 1206 } 1207 }), list( .link = link, .estimated.dispersion = estimated.dispersion, 1208 .type.fitted = type.fitted, 1209 .percentiles = percentiles, 1210 .bred = bred, 1211 .imethod = imethod, .imu = imu, .earg = earg))), 1212 linkinv = eval(substitute(function(eta, extra = NULL) { 1213 mupo <- eta2theta(eta, link = .link , earg = .earg ) 1214 1215 type.fitted <- 1216 if (length(extra$type.fitted)) { 1217 extra$type.fitted 1218 } else { 1219 warning("cannot find 'type.fitted'. Returning 'mean'.") 1220 "mean" 1221 } 1222 1223 type.fitted <- match.arg(type.fitted, 1224 c("mean", "quantiles"))[1] 1225 1226 if (type.fitted == "mean") { 1227 return(label.cols.y(mupo, colnames.y = extra$colnames.y, 1228 NOS = NOS)) 1229 } 1230 1231 1232 percvec <- extra$percentiles 1233 lenperc <- length(percvec) 1234 NOS <- NCOL(eta) / c(M1 = 1) 1235 jvec <- lenperc * (0:(NOS - 1)) 1236 ans <- matrix(0, NROW(eta), lenperc * NOS) 1237 for (kay in 1:lenperc) 1238 ans[, jvec + kay] <- 1239 qpois(0.01 * percvec[kay], lambda = mupo) 1240 1241 rownames(ans) <- rownames(eta) 1242 label.cols.y(ans, colnames.y = extra$colnames.y, 1243 NOS = NOS, percentiles = percvec, 1244 one.on.one = FALSE) 1245 }, list( .link = link, .earg = earg))), 1246 1247 last = eval(substitute(expression({ 1248 if (exists("CQO.FastAlgorithm", envir = VGAMenv)) 1249 rm("CQO.FastAlgorithm", envir = VGAMenv) 1250 dpar <- .dispersion 1251 if (!dpar) { 1252 temp87 <- (y-mu)^2 * 1253 wz / (dtheta.deta(mu, link = .link , earg = .earg )^2) # w cancel 1254 if (M > 1 && ! .onedpar ) { 1255 dpar <- rep_len(NA_real_, M) 1256 temp87 <- cbind(temp87) 1257 nrow.mu <- if (is.matrix(mu)) nrow(mu) else length(mu) 1258 for (ii in 1:M) 1259 dpar[ii] <- sum(temp87[, ii]) / (nrow.mu - ncol(x)) 1260 if (is.matrix(y) && length(dimnames(y)[[2]]) == length(dpar)) 1261 names(dpar) <- dimnames(y)[[2]] 1262 } else { 1263 dpar <- sum(temp87) / (length(mu) - ncol(x)) 1264 } 1265 } 1266 misc$dispersion <- dpar 1267 misc$default.dispersion <- 1 1268 misc$estimated.dispersion <- .estimated.dispersion 1269 1270 misc$imethod <- .imethod 1271 1272 1273 misc$link <- rep_len( .link , M) 1274 names(misc$link) <- if (M > 1) dn2 else new.name # Was old.name=="mu" 1275 1276 misc$earg <- vector("list", M) 1277 names(misc$earg) <- names(misc$link) 1278 for (ii in 1:M) 1279 misc$earg[[ii]] <- .earg 1280 1281 }), list( .dispersion = dispersion, .imethod = imethod, 1282 .estimated.dispersion = estimated.dispersion, 1283 .bred = bred, 1284 .onedpar = onedpar, .link = link, .earg = earg))), 1285 1286 linkfun = eval(substitute( function(mu, extra = NULL) { 1287 theta2eta(mu, link = .link , earg = .earg ) 1288 }, list( .link = link, .earg = earg))), 1289 1290 loglikelihood = eval(substitute( 1291 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 1292 summation = TRUE) { 1293 mupo <- eta2theta(eta, link = .link , earg = .earg ) 1294 if (residuals) { 1295 c(w) * (y / mupo - 1) 1296 } else { 1297 ll.elts <- c(w) * dpois(x = y, lambda = mupo, log = TRUE) 1298 if (summation) { 1299 sum(ll.elts) 1300 } else { 1301 ll.elts 1302 } 1303 } 1304 }, list( .link = link, .earg = earg ))), 1305 vfamily = c("poissonff", "VGAMglm"), # For "stdres" 1306 validparams = eval(substitute(function(eta, y, extra = NULL) { 1307 mupo <- eta2theta(eta, link = .link , earg = .earg ) 1308 okay1 <- all(is.finite(mupo)) && all(0 < mupo) 1309 okay1 1310 }, list( .link = link, .earg = earg ))), 1311 1312 1313 1314 hadof = eval(substitute( 1315 function(eta, extra = list(), deriv = 1, 1316 linpred.index = 1, 1317 w = 1, dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2), 1318 ...) { 1319 mupo <- eta2theta(eta, link = .link , earg = .earg ) 1320 1321 ans <- c(w) * 1322 switch(as.character(deriv), 1323 "0" = 1 / mupo, 1324 "1" = -1 / mupo^2, 1325 "2" = 2 / mupo^3, 1326 "3" = -6 / mupo^4, 1327 stop("argument 'deriv' must be 0, 1, 2 or 3")) 1328 if (deriv == 0) ans else retain.col(ans, linpred.index) # Coz M1 = 1 1329 }, list( .link = link, .earg = earg ))), 1330 1331 1332 1333 1334 1335 simslot = 1336 function(object, nsim) { 1337 1338 1339 pwts <- if (length(pwts <- object@prior.weights) > 0) 1340 pwts else weights(object, type = "prior") 1341 if (any(pwts != 1)) 1342 warning("ignoring prior weights") 1343 ftd <- fitted(object) 1344 rpois(nsim * length(ftd), ftd) 1345 }, 1346 1347 1348 1349 deriv = eval(substitute(expression({ 1350 mupo <- eta2theta(eta, link = .link , earg = .earg ) 1351 yBRED <- if ( .bred ) { 1352 Hvector <- hatvaluesbasic(X.vlm = X.vlm.save, 1353 diagWm = c(t(c(w) * mupo))) # Handles M>1 1354 1355 1356 varY <- mupo # Is a matrix if M>1. 1357 d1.BRED <- dtheta.deta(mupo, .link , earg = .earg ) 1358 d2.BRED <- d2theta.deta2(mupo, .link , earg = .earg ) 1359 y + matrix(Hvector, n, M, byrow = TRUE) * 1360 varY * d2.BRED / (2 * d1.BRED^2) 1361 } else { 1362 y 1363 } 1364 1365 1366 answer <- if ( .link == "loglink" && (any(mupo < .Machine$double.eps))) { 1367 c(w) * (yBRED - mupo) 1368 } else { 1369 lambda <- mupo 1370 dl.dlambda <- (yBRED - lambda) / lambda 1371 dlambda.deta <- dtheta.deta(theta = lambda, 1372 link = .link , earg = .earg ) 1373 c(w) * dl.dlambda * dlambda.deta 1374 } 1375 1376 answer 1377 }), list( .link = link, .earg = earg, .bred = bred))), 1378 1379 weight = eval(substitute(expression({ 1380 if ( .link == "loglink" && (any(mupo < .Machine$double.eps))) { 1381 tmp600 <- mupo 1382 tmp600[tmp600 < .Machine$double.eps] <- .Machine$double.eps 1383 c(w) * tmp600 1384 } else { 1385 ned2l.dlambda2 <- 1 / lambda 1386 ned2lambda.deta2 <- d2theta.deta2(theta = lambda, 1387 link = .link , earg = .earg ) 1388 c(w) * dlambda.deta^2 * ned2l.dlambda2 1389 } 1390 }), list( .link = link, .earg = earg)))) 1391} # poissonff() 1392 1393 1394 1395 1396if (FALSE) 1397 quasibinomialff <- 1398 function( 1399 link = "logitlink", 1400 multiple.responses = FALSE, onedpar = !multiple.responses, 1401 parallel = FALSE, zero = NULL) { 1402 1403 1404 link <- as.list(substitute(link)) 1405 earg <- link2list(link) 1406 link <- attr(earg, "function.name") 1407 1408 dispersion <- 0 # Estimated; this is the only difference w. binomialff() 1409 ans <- binomialff(link = earg, earg.link = TRUE, 1410 dispersion = dispersion, 1411 multiple.responses = multiple.responses, 1412 onedpar = onedpar, 1413 parallel = parallel, zero = zero) 1414 ans@vfamily <- "quasibinomialff" 1415 ans@infos <- eval(substitute(function(...) { 1416 list(M1 = 1, 1417 Q1 = 1, 1418 multipleResponses = .multiple.responses , 1419 parameters.names = c("prob"), 1420 quasi.type = TRUE, 1421 zero = .zero ) 1422 }, list( .zero = zero, 1423 .multiple.responses = multiple.responses ))) 1424 1425 ans 1426} # quasibinomialff 1427 1428 1429 1430 1431 1432if (FALSE) 1433 quasipoissonff <- function(link = "loglink", onedpar = FALSE, 1434 parallel = FALSE, zero = NULL) { 1435 1436 link <- as.list(substitute(link)) 1437 earg <- link2list(link) 1438 link <- attr(earg, "function.name") 1439 1440 1441 1442 dispersion <- 0 # Ested; this is the only difference with poissonff() 1443 ans <- poissonff(link = earg, earg.link = TRUE, 1444 dispersion = dispersion, onedpar = onedpar, 1445 parallel = parallel, zero = zero) 1446 ans@vfamily <- "quasipoissonff" 1447 ans@infos <- eval(substitute(function(...) { 1448 list(M1 = 1, 1449 Q1 = 1, 1450 multipleResponses = TRUE, 1451 parameters.names = c("lambda"), 1452 quasi.type = TRUE, 1453 zero = .zero ) 1454 }, list( .zero = zero ))) 1455 1456 ans 1457} # quasipoissonff 1458 1459 1460 1461 1462 1463 1464 double.exppoisson <- 1465 function(lmean = "loglink", 1466 ldispersion = "logitlink", 1467 idispersion = 0.8, 1468 zero = NULL) { 1469 1470 if (!is.Numeric(idispersion, positive = TRUE)) 1471 stop("bad input for 'idispersion'") 1472 1473 1474 lmean <- as.list(substitute(lmean)) 1475 emean <- link2list(lmean) 1476 lmean <- attr(emean, "function.name") 1477 1478 ldisp <- as.list(substitute(ldispersion)) 1479 edisp <- link2list(ldisp) 1480 ldisp <- attr(edisp, "function.name") 1481 1482 idisp <- idispersion 1483 1484 1485 new("vglmff", 1486 blurb = c("Double exponential Poisson distribution\n\n", 1487 "Link: ", 1488 namesof("mean", lmean, earg = emean), ", ", 1489 namesof("dispersion", ldisp, earg = edisp), "\n", 1490 "Mean: ", "mean\n", 1491 "Variance: mean / dispersion"), 1492 constraints = eval(substitute(expression({ 1493 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1494 predictors.names = predictors.names, 1495 M1 = 2) 1496 }), list( .zero = zero ))), 1497 1498 infos = eval(substitute(function(...) { 1499 list(M1 = 2, 1500 parameters.names = c("mean", "dispersion"), 1501 lmean = .lmean , 1502 ldispersion = .ldispersion , 1503 zero = .zero ) 1504 }, list( .lmean = lmean, 1505 .ldispersion = ldispersion, 1506 .zero = zero ))), 1507 1508 1509 initialize = eval(substitute(expression({ 1510 1511 w.y.check(w = w, y = y, 1512 Is.nonnegative.y = TRUE, 1513 ncol.w.max = 1, 1514 ncol.y.max = 1) 1515 1516 1517 M <- if (is.matrix(y)) ncol(y) else 1 1518 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 1519 dn2 <- if (length(dn2)) { 1520 paste("E[", dn2, "]", sep = "") 1521 } else { 1522 "mu" 1523 } 1524 predictors.names <- 1525 c(namesof(dn2, .lmean , earg = .emean, short = TRUE), 1526 namesof("dispersion", .ldisp , earg = .edisp, short = TRUE)) 1527 1528 init.mu <- pmax(y, 1/8) 1529 tmp2 <- rep_len( .idisp , n) 1530 1531 if (!length(etastart)) 1532 etastart <- 1533 cbind(theta2eta(init.mu, .lmean , earg = .emean ), 1534 theta2eta(tmp2, .ldisp , earg = .edisp )) 1535 }), list( .lmean = lmean, .emean = emean, 1536 .ldisp = ldisp, .edisp = edisp, 1537 .idisp = idisp ))), 1538 linkinv = eval(substitute(function(eta, extra = NULL) { 1539 eta2theta(eta[, 1], link = .lmean, earg = .emean) 1540 }, list( .lmean = lmean, .emean = emean, 1541 .ldisp = ldisp, .edisp = edisp ))), 1542 last = eval(substitute(expression({ 1543 misc$link <- c(mean = .lmean , dispersion = .ldisp ) 1544 1545 misc$earg <- list(mean = .emean , dispersion = .edisp ) 1546 1547 misc$expected <- TRUE 1548 }), list( .lmean = lmean, .emean = emean, 1549 .ldisp = ldisp, .edisp = edisp ))), 1550 loglikelihood = eval(substitute( 1551 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 1552 summation = TRUE) { 1553 lambda <- eta2theta(eta[, 1], .lmean , earg = .emean ) 1554 Disper <- eta2theta(eta[, 2], .ldisp , earg = .edisp ) 1555 if (residuals) { 1556 stop("loglikelihood residuals not implemented yet") 1557 } else { 1558 ll.elts <- c(w) * (0.5 * log(Disper) + 1559 Disper*(y-lambda) + Disper*y*log(lambda)) 1560 if (summation) { 1561 sum(ll.elts) 1562 } else { 1563 ll.elts 1564 } 1565 } 1566 }, list( .lmean = lmean, .emean = emean, 1567 .ldisp = ldisp, .edisp = edisp ))), 1568 vfamily = "double.exppoisson", 1569 validparams = eval(substitute(function(eta, y, extra = NULL) { 1570 lambda <- eta2theta(eta[, 1], .lmean , earg = .emean ) 1571 Disper <- eta2theta(eta[, 2], .ldisp , earg = .edisp ) 1572 okay1 <- all(is.finite(lambda)) && all(0 < lambda) && 1573 all(is.finite(Disper)) && all(0 < Disper & Disper < 1) 1574 okay1 1575 }, list( .lmean = lmean, .emean = emean, 1576 .ldisp = ldisp, .edisp = edisp ))), 1577 1578 1579 deriv = eval(substitute(expression({ 1580 lambda <- eta2theta(eta[, 1], .lmean , earg = .emean ) 1581 Disper <- eta2theta(eta[, 2], .ldisp , earg = .edisp ) 1582 1583 dl.dlambda <- Disper * (y / lambda - 1) 1584 dl.dDisper <- y * log(lambda) + y - lambda + 0.5 / Disper 1585 1586 dlambda.deta <- dtheta.deta(theta = lambda, .lmean, earg = .emean) 1587 dDisper.deta <- dtheta.deta(theta = Disper, .ldisp, earg = .edisp) 1588 1589 c(w) * cbind(dl.dlambda * dlambda.deta, 1590 dl.dDisper * dDisper.deta) 1591 }), list( .lmean = lmean, .emean = emean, 1592 .ldisp = ldisp, .edisp = edisp ))), 1593 weight = eval(substitute(expression({ 1594 wz <- matrix(NA_real_, nrow = n, ncol = 2) # diagonal 1595 usethis.lambda <- pmax(lambda, .Machine$double.eps / 10000) 1596 wz[, iam(1, 1, M)] <- (Disper / usethis.lambda) * dlambda.deta^2 1597 wz[, iam(2, 2, M)] <- (0.5 / Disper^2) * dDisper.deta^2 1598 c(w) * wz 1599 }), list( .lmean = lmean, .emean = emean, 1600 .ldisp = ldisp, 1601 .edisp = edisp )))) 1602} # double.exppoisson 1603 1604 1605 1606 double.expbinomial <- 1607 function(lmean = "logitlink", ldispersion = "logitlink", 1608 idispersion = 0.25, zero = "dispersion") { 1609 1610 lmean <- as.list(substitute(lmean)) 1611 emean <- link2list(lmean) 1612 lmean <- attr(emean, "function.name") 1613 1614 ldisp <- as.list(substitute(ldispersion)) 1615 edisp <- link2list(ldisp) 1616 ldisp <- attr(edisp, "function.name") 1617 idisp <- idispersion 1618 1619 1620 if (!is.Numeric(idispersion, positive = TRUE)) 1621 stop("bad input for 'idispersion'") 1622 1623 1624 new("vglmff", 1625 blurb = c("Double Exponential Binomial distribution\n\n", 1626 "Link: ", 1627 namesof("mean", lmean, earg = emean), ", ", 1628 namesof("dispersion", ldisp, earg = edisp), "\n", 1629 "Mean: ", "mean\n"), 1630 constraints = eval(substitute(expression({ 1631 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1632 predictors.names = predictors.names, 1633 M1 = 2) 1634 }), list( .zero = zero ))), 1635 1636 1637 infos = eval(substitute(function(...) { 1638 list(M1 = 2, 1639 Q1 = NA, 1640 parameters.names = c("mean", "dispersion"), 1641 lmean = .lmean , 1642 ldisp = .ldisp , 1643 multipleResponses = FALSE, 1644 zero = .zero ) 1645 }, list( .lmean = lmean, 1646 .zero = zero, 1647 .ldisp = ldisp ))), 1648 1649 initialize = eval(substitute(expression({ 1650 if (!all(w == 1)) 1651 extra$orig.w <- w 1652 1653 1654 if (NCOL(w) != 1) 1655 stop("'weights' must be a vector or a one-column matrix") 1656 1657 NCOL <- function (x) 1658 if (is.array(x) && length(dim(x)) > 1 || 1659 is.data.frame(x)) ncol(x) else as.integer(1) 1660 1661 if (NCOL(y) == 1) { 1662 1663 1664 if (is.factor(y)) y <- (y != levels(y)[1]) 1665 nvec <- rep_len(1, n) 1666 y[w == 0] <- 0 1667 if (!all(y == 0 | y == 1)) 1668 stop("response values 'y' must be 0 or 1") 1669 init.mu <- (0.5 + w * y) / (1 + w) 1670 1671 1672 no.successes <- y 1673 if (min(y) < 0) 1674 stop("Negative data not allowed!") 1675 if (any(abs(no.successes - round(no.successes)) > 1.0e-8)) 1676 stop("Number of successes must be integer-valued") 1677 } else if (NCOL(y) == 2) { 1678 if (min(y) < 0) 1679 stop("Negative data not allowed!") 1680 if (any(abs(y - round(y)) > 1.0e-8)) 1681 stop("Count data must be integer-valued") 1682 y <- round(y) 1683 nvec <- y[, 1] + y[, 2] 1684 y <- ifelse(nvec > 0, y[, 1] / nvec, 0) 1685 w <- w * nvec 1686 init.mu <- (0.5 + nvec * y) / (1 + nvec) 1687 } else 1688 stop("for the double.expbinomial family, response 'y' must be", 1689 " a vector of 0 and 1's\n", 1690 "or a factor (first level = fail, ", 1691 "other levels = success),\n", 1692 "or a 2-column matrix where col 1 is the no. of ", 1693 "successes and col 2 is the no. of failures") 1694 1695 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 1696 dn2 <- if (length(dn2)) paste("E[", dn2, "]", sep = "") else "mu" 1697 1698 predictors.names <- 1699 c(namesof(dn2, .lmean , earg = .emean , short = TRUE), 1700 namesof("dispersion", .ldisp , earg = .edisp , short = TRUE)) 1701 1702 tmp2 <- rep_len( .idisp , n) 1703 1704 if (!length(etastart)) 1705 etastart <- cbind(theta2eta(init.mu, .lmean, earg = .emean), 1706 theta2eta(tmp2, .ldisp, earg = .edisp)) 1707 }), list( .lmean = lmean, .emean = emean, 1708 .ldisp = ldisp, .edisp = edisp, 1709 .idisp = idisp ))), 1710 linkinv = eval(substitute(function(eta, extra = NULL) { 1711 eta2theta(eta[, 1], link = .lmean , earg = .emean ) 1712 }, list( .lmean = lmean, .emean = emean, 1713 .ldisp = ldisp, .edisp = edisp ))), 1714 last = eval(substitute(expression({ 1715 misc$link <- c("mean" = .lmean, "dispersion" = .ldisp) 1716 1717 misc$earg <- list( mean = .emean, dispersion = .edisp) 1718 1719 misc$expected <- TRUE 1720 }), list( .lmean = lmean, .emean = emean, 1721 .ldisp = ldisp, .edisp = edisp ))), 1722 loglikelihood = eval(substitute( 1723 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 1724 summation = TRUE) { 1725 prob <- eta2theta(eta[, 1], link = .lmean, earg = .emean) 1726 Disper <- eta2theta(eta[, 2], link = .ldisp, earg = .edisp) 1727 if (residuals) { 1728 stop("loglikelihood residuals not implemented yet") 1729 } else { 1730 1731 1732 1733 temp1 <- y * log(ifelse(y > 0, y, 1)) # y*log(y) 1734 temp2 <- (1.0-y)* log1p(ifelse(y < 1, -y, 0)) # (1-y)*log(1-y) 1735 1736 1737 ll.elts <- 1738 (0.5 * log(Disper) + w * (y * Disper * log(prob) + 1739 (1-y) * Disper * log1p(-prob) + 1740 temp1 * (1-Disper) + temp2 * (1 - Disper))) 1741 if (summation) { 1742 sum(ll.elts) 1743 } else { 1744 ll.elts 1745 } 1746 } 1747 }, list( .lmean = lmean, .emean = emean, 1748 .ldisp = ldisp, .edisp = edisp ))), 1749 vfamily = "double.expbinomial", 1750 validparams = eval(substitute(function(eta, y, extra = NULL) { 1751 prob <- eta2theta(eta[, 1], link = .lmean , earg = .emean ) 1752 Disper <- eta2theta(eta[, 2], link = .ldisp , earg = .edisp ) 1753 okay1 <- all(is.finite(prob )) && all(0 < prob & prob < 1) && 1754 all(is.finite(Disper)) && all(0 < Disper & Disper < 1) 1755 okay1 1756 }, list( .lmean = lmean, .emean = emean, 1757 .ldisp = ldisp, .edisp = edisp ))), 1758 deriv = eval(substitute(expression({ 1759 prob <- eta2theta(eta[, 1], link = .lmean , earg = .emean ) 1760 Disper <- eta2theta(eta[, 2], link = .ldisp , earg = .edisp ) 1761 temp1 <- y * log(ifelse(y > 0, y, 1)) # y*log(y) 1762 temp2 <- (1.0-y) * log1p(ifelse(y < 1, -y, 0)) # (1-y)*log(1-y) 1763 temp3 <- prob * (1.0-prob) 1764 temp3 <- pmax(temp3, .Machine$double.eps * 10000) 1765 1766 dl.dprob <- w * Disper * (y - prob) / temp3 1767 dl.dDisper <- 0.5 / Disper + w * (y * log(prob) + 1768 (1-y)*log1p(-prob) - temp1 - temp2) 1769 1770 dprob.deta <- dtheta.deta(theta = prob, .lmean , earg = .emean ) 1771 dDisper.deta <- dtheta.deta(theta = Disper, .ldisp , earg = .edisp ) 1772 1773 cbind(dl.dprob * dprob.deta, 1774 dl.dDisper * dDisper.deta) 1775 }), list( .lmean = lmean, .emean = emean, 1776 .ldisp = ldisp, .edisp = edisp ))), 1777 weight = eval(substitute(expression({ 1778 wz <- matrix(NA_real_, nrow = n, ncol = 2) # diagonal 1779 wz[, iam(1, 1, M)] <- w * (Disper / temp3) * dprob.deta^2 1780 wz[, iam(2, 2, M)] <- (0.5 / Disper^2) * dDisper.deta^2 1781 wz 1782 }), list( .lmean = lmean, .emean = emean, 1783 .ldisp = ldisp, .edisp = edisp )))) 1784} # double.expbinomial 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 augbinomial <- function(link = "logitlink", multiple.responses = FALSE, 1796 parallel = TRUE) { 1797 1798 if (!is.logical(parallel) || 1799 length(parallel) != 1 || 1800 !parallel) 1801 warning("Argument 'parallel' should be assigned 'TRUE' only") 1802 1803 link <- as.list(substitute(link)) 1804 earg <- link2list(link) 1805 link <- attr(earg, "function.name") 1806 1807 1808 new("vglmff", 1809 blurb = if (multiple.responses) 1810 c("Augmented multivariate binomial model\n\n", 1811 "Link: ", 1812 namesof("mu.1[,j]", link, earg = earg), ", ", 1813 namesof("mu.2[,j]", link, earg = earg), 1814 "\n", 1815 "Variance: mu[,j]*(1-mu[,j])") else 1816 c("Augmented binomial model\n\n", 1817 "Link: ", 1818 namesof("mu.1[,j]", link, earg = earg), ", ", 1819 namesof("mu.2[,j]", link, earg = earg), 1820 "\n", 1821 "Variance: mu*(1-mu)"), 1822 deviance = function(mu, y, w, residuals = FALSE, eta, extra = NULL, 1823 summation = TRUE) { 1824 Deviance.categorical.data.vgam(mu = cbind(mu, 1-mu), y=cbind(y, 1-y), 1825 w = w, residuals = residuals, 1826 eta = eta, extra = extra, 1827 summation = summation) 1828 }, 1829 infos = eval(substitute(function(...) { 1830 list(M1 = 2, 1831 parameters.names = c("mu.1[,j]", "mu.2[,j]"), 1832 parallel = .parallel) 1833 }, list( .parallel = parallel ))), 1834 initialize = eval(substitute(expression({ 1835 1836 M1 = 2 1837 1838 if ( .multiple.responses ) { 1839 y <- as.matrix(y) 1840 M <- M1 * ncol(y) 1841 if (!all(y == 0 | y == 1)) 1842 stop("response must contain 0's and 1's only") 1843 dn2 <- if (is.matrix(y)) dimnames(y)[[2]] else NULL 1844 dn2 <- if (length(dn2)) { 1845 paste("E[", dn2, "]", sep = "") 1846 } else { 1847 param.names("mu", M, skip1 = TRUE) 1848 } 1849 predictors.names <- 1850 c(namesof(if (M > 1) dn2 else 1851 "mu.1", .link , earg = .earg , short = TRUE), 1852 namesof(if (M > 1) dn2 else 1853 "mu.2", .link , earg = .earg , short = TRUE)) 1854 NOS <- M / M1 1855 predictors.names <- 1856 predictors.names[interleave.VGAM(M1 * NOS, M1 = M1)] 1857 1858 1859 if (!length(mustart) && !length(etastart)) 1860 mustart <- (0.5 + w * y) / (1 + w) 1861 } else { 1862 1863 dn2 <- c("mu1.", "mu2.") 1864 M <- M1 1865 1866 1867 1868 if (!all(w == 1)) 1869 extra$orig.w <- w 1870 1871 1872 NCOL <- function (x) if (is.array(x) && length(dim(x)) > 1 || 1873 is.data.frame(x)) ncol(x) else as.integer(1) 1874 if (NCOL(y) == 1) { 1875 if (is.factor(y)) y = (y != levels(y)[1]) 1876 nvec <- rep_len(1, n) 1877 y[w == 0] <- 0 1878 if (!all(y == 0 | y == 1)) 1879 stop("response values 'y' must be 0 or 1") 1880 if (!length(mustart) && !length(etastart)) 1881 mustart <- (0.5 + w * y) / (1 + w) 1882 1883 1884 no.successes <- y 1885 if (min(y) < 0) 1886 stop("Negative data not allowed!") 1887 if (any(abs(no.successes - round(no.successes)) > 1.0e-8)) 1888 stop("Number of successes must be integer-valued") 1889 } else if (NCOL(y) == 2) { 1890 if (min(y) < 0) 1891 stop("Negative data not allowed!") 1892 if (any(abs(y - round(y)) > 1.0e-8)) 1893 stop("Count data must be integer-valued") 1894 y <- round(y) 1895 nvec <- y[, 1] + y[, 2] 1896 y <- ifelse(nvec > 0, y[, 1] / nvec, 0) 1897 w <- w * nvec 1898 if (!length(mustart) && !length(etastart)) 1899 mustart <- (0.5 + nvec * y) / (1 + nvec) 1900 } else { 1901 stop("for the binomialff family, response 'y' must be a ", 1902 "vector of 0 and 1's\n", 1903 "or a factor (first level = fail, ", 1904 "other levels = success),\n", 1905 "or a 2-column matrix where col 1 is the no. of ", 1906 "successes and col 2 is the no. of failures") 1907 } 1908 predictors.names <- 1909 c(namesof("mu.1", .link , earg = .earg , short = TRUE), 1910 namesof("mu.2", .link , earg = .earg , short = TRUE)) 1911 } 1912 }), list( .link = link, 1913 .multiple.responses = multiple.responses, .earg = earg))), 1914 linkinv = eval(substitute(function(eta, extra = NULL) { 1915 Mdiv2 <- ncol(eta) / 2 1916 index1 <- 2*(1:Mdiv2) - 1 1917 mu <- eta2theta(eta[, index1], link = .link , earg = .earg ) 1918 mu 1919 }, list( .link = link, .earg = earg ))), 1920 last = eval(substitute(expression({ 1921 misc$link <- rep_len( .link , M) 1922 names(misc$link) <- if (M > 1) dn2 else "mu" 1923 1924 misc$earg <- vector("list", M) 1925 names(misc$earg) <- names(misc$link) 1926 for (ii in 1:M) 1927 misc$earg[[ii]] <- .earg 1928 1929 misc$parallel <- .parallel 1930 misc$expected <- TRUE 1931 misc$multiple.responses <- .multiple.responses 1932 }), list( .link = link, 1933 .multiple.responses = multiple.responses, .earg = earg, 1934 .parallel = parallel ))), 1935 linkfun = eval(substitute(function(mu, extra = NULL) { 1936 usualanswer <- theta2eta(mu, .link , earg = .earg ) 1937 kronecker(usualanswer, matrix(1, 1, 2)) 1938 }, list( .link = link, .earg = earg))), 1939 loglikelihood = 1940 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 1941 summation = TRUE) { 1942 if (residuals) { 1943 c(w) * (y / mu - (1-y) / (1-mu)) 1944 } else { 1945 ycounts <- if (is.numeric(extra$orig.w)) y * w / extra$orig.w else 1946 y * c(w) # Convert proportions to counts 1947 nvec <- if (is.numeric(extra$orig.w)) round(w / extra$orig.w) else 1948 round(w) 1949 1950 smallno <- 1.0e6 * .Machine$double.eps 1951 smallno <- sqrt(.Machine$double.eps) 1952 if (max(abs(ycounts - round(ycounts))) > smallno) 1953 warning("converting 'ycounts' to integer in @loglikelihood") 1954 ycounts <- round(ycounts) 1955 1956 ll.elts <- 1957 (if (is.numeric(extra$orig.w)) extra$orig.w else 1) * 1958 dbinom(x = ycounts, size = nvec, prob = mu, log = TRUE) 1959 if (summation) { 1960 sum(ll.elts) 1961 } else { 1962 ll.elts 1963 } 1964 } 1965 }, 1966 vfamily = c("augbinomial", "VGAMcategorical"), 1967 validparams = eval(substitute(function(eta, y, extra = NULL) { 1968 Mdiv2 <- ncol(eta) / 2 1969 index1 <- 2*(1:Mdiv2) - 1 1970 mu <- eta2theta(eta[, index1], link = .link , earg = .earg ) 1971 okay1 <- all(is.finite(mu)) && all(0 < mu & mu < 1) 1972 okay1 1973 }, list( .link = link, .earg = earg))), 1974 deriv = eval(substitute(expression({ 1975 M1 <- 2 1976 Mdiv2 <- M / 2 1977 NOS <- M / M1 1978 1979 Konst1 <- 1 # Works with this 1980 deriv1 <- Konst1 * w * 1981 if ( .link == "logitlink") { 1982 y * (1 - mu) 1983 } else { 1984 stop("this is not programmed in yet") 1985 dtheta.deta(mu, link = .link , earg = .earg ) * 1986 (y / mu - 1.0) / (1.0 - mu) 1987 } 1988 deriv2 = Konst1 * w * 1989 if ( .link == "logitlink") { 1990 -(1 - y) * mu 1991 } else { 1992 stop("this is not programmed in yet") 1993 dtheta.deta(mu, link = .link , earg = .earg ) * 1994 (y / mu - 1.0) / (1.0 - mu) 1995 } 1996 1997 myderiv <- (cbind(deriv1, 1998 deriv2))[, interleave.VGAM(M1 * NOS, M1 = M1)] 1999 myderiv 2000 }), list( .link = link, .earg = earg))), 2001 weight = eval(substitute(expression({ 2002 tmp100 <- mu * (1.0 - mu) 2003 2004 tmp200 <- if ( .link == "logitlink") { 2005 cbind(w * tmp100) 2006 } else { 2007 cbind(w * dtheta.deta(mu, link = .link , earg = .earg )^2 / tmp100) 2008 } 2009 2010 wk.wt1 <- (Konst1^2) * tmp200 * (1 - mu) 2011 wk.wt2 <- (Konst1^2) * tmp200 * mu 2012 2013 2014 2015 2016 my.wk.wt <- cbind(wk.wt1, wk.wt2) 2017 my.wk.wt <- my.wk.wt[, interleave.VGAM(M1 * NOS, M1 = M1)] 2018 my.wk.wt 2019 }), list( .link = link, .earg = earg)))) 2020} 2021 2022 2023 2024 2025 2026