1# These functions are 2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland. 3# All rights reserved. 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28dlms.bcn <- function(x, lambda = 1, mu = 0, sigma = 1, 29 tol0 = 0.001, log = FALSE) { 30 if (!is.logical(log.arg <- log) || length(log) != 1) 31 stop("bad input for argument 'log'") 32 rm(log) 33 34 zedd <- ((x/mu)^lambda - 1) / (lambda * sigma) 35 log.dz.dy <- (lambda - 1) * log(x/mu) - log(mu * sigma) 36 37 is.eff.0 <- abs(lambda) < tol0 38 if (any(is.eff.0)) { 39 zedd[is.eff.0] <- log(x[is.eff.0] / mu[is.eff.0]) / sigma[is.eff.0] 40 log.dz.dy[is.eff.0] <- -log(x[is.eff.0] * sigma[is.eff.0]) 41 } 42 logden <- dnorm(zedd, log = TRUE) + log.dz.dy 43 if (log.arg) logden else exp(logden) 44} 45 46 47 48qlms.bcn <- function(p, lambda = 1, mu = 0, sigma = 1) { 49 50 answer <- mu * (1 + lambda * sigma * qnorm(p))^(1/lambda) 51 answer 52} 53 54 55 56 57 58 59 60lms.bcn.control <- 61lms.bcg.control <- 62lms.yjn.control <- function(trace = TRUE, ...) 63 list(trace = trace) 64 65 66 67 lms.bcn <- function(percentiles = c(25, 50, 75), 68 zero = c("lambda", "sigma"), 69 llambda = "identitylink", 70 lmu = "identitylink", 71 lsigma = "loglink", 72 idf.mu = 4, 73 idf.sigma = 2, 74 ilambda = 1, 75 isigma = NULL, 76 tol0 = 0.001) { 77 llambda <- as.list(substitute(llambda)) 78 elambda <- link2list(llambda) 79 llambda <- attr(elambda, "function.name") 80 81 82 lmu <- as.list(substitute(lmu)) 83 emu <- link2list(lmu) 84 lmu <- attr(emu, "function.name") 85 86 lsigma <- as.list(substitute(lsigma)) 87 esigma <- link2list(lsigma) 88 lsigma <- attr(esigma, "function.name") 89 90 91 if (!is.Numeric(tol0, positive = TRUE, length.arg = 1)) 92 stop("bad input for argument 'tol0'") 93 94 if (!is.Numeric(ilambda)) 95 stop("bad input for argument 'ilambda'") 96 if (length(isigma) && 97 !is.Numeric(isigma, positive = TRUE)) 98 stop("bad input for argument 'isigma'") 99 100 101 102 new("vglmff", 103 blurb = c("LMS ", 104 "quantile", 105 " regression (Box-Cox transformation to normality)\n", 106 "Links: ", 107 namesof("lambda", link = llambda, earg = elambda), ", ", 108 namesof("mu", link = lmu, earg = emu), ", ", 109 namesof("sigma", link = lsigma, earg = esigma)), 110 constraints = eval(substitute(expression({ 111 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 112 predictors.names = predictors.names, 113 M1 = 3) 114 }), list( .zero = zero))), 115 116 infos = eval(substitute(function(...) { 117 list(M1 = 3, 118 Q1 = 1, 119 expected = TRUE, 120 multipleResponses = FALSE, 121 parameters.names = c("lambda", "mu", "sigma"), 122 llambda = .llambda , 123 lmu = .lmu , 124 lsigma = .lsigma , 125 percentiles = .percentiles , # For the original fit only 126 true.mu = FALSE, # quantiles 127 zero = .zero ) 128 }, list( .zero = zero, 129 .percentiles = percentiles, 130 .llambda = llambda, .lmu = lmu, .lsigma = lsigma ))), 131 132 initialize = eval(substitute(expression({ 133 134 w.y.check(w = w, y = y, 135 Is.positive.y = TRUE, 136 ncol.w.max = 1, ncol.y.max = 1) 137 138 139 predictors.names <- 140 c(namesof("lambda", .llambda, earg = .elambda, short= TRUE), 141 namesof("mu", .lmu, earg = .emu, short= TRUE), 142 namesof("sigma", .lsigma, earg = .esigma, short= TRUE)) 143 144 extra$percentiles <- .percentiles 145 146 147 if (!length(etastart)) { 148 149 Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 150 y = y, w = w, df = .idf.mu ) 151 fv.init <- c(predict(Fit5, x = x[, min(ncol(x), 2)])$y) 152 153 lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0 154 sigma.init <- if (is.null(.isigma)) { 155 myratio <- ((y/fv.init)^lambda.init - 1) / lambda.init 156 if (is.Numeric( .idf.sigma )) { 157 fit600 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 158 y = myratio^2, 159 w = w, df = .idf.sigma) 160 sqrt(c(abs(predict(fit600, x = x[, min(ncol(x), 2)])$y))) 161 } else { 162 sqrt(var(myratio)) 163 } 164 } else { 165 .isigma 166 } 167 168 etastart <- 169 cbind(theta2eta(lambda.init, .llambda , earg = .elambda ), 170 theta2eta(fv.init, .lmu , earg = .emu ), 171 theta2eta(sigma.init, .lsigma , earg = .esigma )) 172 } 173 }), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 174 .elambda = elambda, .emu = emu, .esigma = esigma, 175 .idf.mu = idf.mu, 176 .idf.sigma = idf.sigma, 177 .ilambda = ilambda, .isigma = isigma, 178 .percentiles = percentiles ))), 179 180 181 linkinv = eval(substitute(function(eta, extra = NULL) { 182 pcent <- extra$percentiles 183 184 eta[, 1] <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 185 eta[, 2] <- eta2theta(eta[, 2], .lmu , earg = .emu ) 186 eta[, 3] <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 187 qtplot.lms.bcn(percentiles = pcent, eta = eta) 188 }, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 189 .elambda = elambda, .emu = emu, .esigma = esigma 190 ))), 191 192 193 last = eval(substitute(expression({ 194 misc$links <- c(lambda = .llambda , mu = .lmu , sigma = .lsigma ) 195 196 misc$earg <- list(lambda = .elambda , mu = .emu , sigma = .esigma ) 197 198 misc$tol0 <- .tol0 199 misc$percentiles <- .percentiles # These are argument values 200 if (control$cdf) { 201 post$cdf <- 202 cdf.lms.bcn(y, 203 eta0 = matrix(c(lambda, mymu, sigma), ncol = 3, 204 dimnames = list(dimnames(x)[[1]], NULL))) 205 } 206 }), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 207 .elambda = elambda, .emu = emu, .esigma = esigma, 208 .percentiles = percentiles, 209 .tol0 = tol0 ))), 210 loglikelihood = eval(substitute( 211 function(mu, y, w, residuals = FALSE, eta, 212 extra = NULL, 213 summation = TRUE) { 214 215 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 216 muvec <- eta2theta(eta[, 2], .lmu , earg = .emu ) 217 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 218 219 220 if (residuals) { 221 stop("loglikelihood residuals not implemented yet") 222 } else { 223 ll.elts <- dlms.bcn(x = y, lambda = lambda, mu = mu, sigma = sigma, 224 tol0 = .tol0 , log = TRUE) 225 if (summation) { 226 sum(ll.elts) 227 } else { 228 ll.elts 229 } 230 } 231 }, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 232 .elambda = elambda, .emu = emu, .esigma = esigma, 233 .tol0 = tol0 ))), 234 vfamily = c("lms.bcn", "lmscreg"), 235 validparams = eval(substitute(function(eta, y, extra = NULL) { 236 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 237 mymu <- eta2theta(eta[, 2], .lmu , earg = .emu ) 238 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 239 okay1 <- all(is.finite(mymu )) && 240 all(is.finite(sigma )) && all(0 < sigma) && 241 all(is.finite(lambda)) 242 okay1 243 }, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 244 .elambda = elambda, .emu = emu, .esigma = esigma, 245 .tol0 = tol0 ))), 246 deriv = eval(substitute(expression({ 247 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 248 mymu <- eta2theta(eta[, 2], .lmu , earg = .emu ) 249 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 250 251 zedd <- ((y / mymu)^lambda - 1) / (lambda * sigma) 252 z2m1 <- zedd * zedd - 1 253 254 dl.dlambda <- zedd * (zedd - log(y/mymu) / sigma) / lambda - 255 z2m1 * log(y/mymu) 256 dl.dmu <- zedd / (mymu * sigma) + z2m1 * lambda / mymu 257 dl.dsigma <- z2m1 / sigma 258 259 dlambda.deta <- dtheta.deta(lambda, .llambda , earg = .elambda ) 260 dmu.deta <- dtheta.deta(mymu, .lmu , earg = .emu ) 261 dsigma.deta <- dtheta.deta(sigma, .lsigma , earg = .esigma ) 262 263 c(w) * cbind(dl.dlambda * dlambda.deta, 264 dl.dmu * dmu.deta, 265 dl.dsigma * dsigma.deta) 266 }), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 267 .elambda = elambda, .emu = emu, .esigma = esigma ))), 268 weight = eval(substitute(expression({ 269 wz <- matrix(NA_real_, n, 6) 270 wz[,iam(1, 1, M)] <- (7 * sigma^2 / 4) * dlambda.deta^2 271 wz[,iam(2, 2, M)] <- (1 + 2*(lambda*sigma)^2)/(mymu*sigma)^2 * 272 dmu.deta^2 273 wz[,iam(3, 3, M)] <- (2 / sigma^2) * dsigma.deta^2 274 wz[,iam(1, 2, M)] <- (-1 / (2 * mymu)) * dlambda.deta * dmu.deta 275 wz[,iam(1, 3, M)] <- (lambda * sigma) * dlambda.deta * dsigma.deta 276 wz[,iam(2, 3, M)] <- (2*lambda/(mymu * sigma)) * 277 dmu.deta * dsigma.deta 278 c(w) * wz 279 }), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 280 .elambda = elambda, .emu = emu, .esigma = esigma )))) 281} # End of lms.bcn 282 283 284 285 286 287 288 lms.bcg <- function(percentiles = c(25, 50, 75), 289 zero = c("lambda", "sigma"), 290 llambda = "identitylink", 291 lmu = "identitylink", 292 lsigma = "loglink", 293 idf.mu = 4, 294 idf.sigma = 2, 295 ilambda = 1, 296 isigma = NULL) { 297 llambda <- as.list(substitute(llambda)) 298 elambda <- link2list(llambda) 299 llambda <- attr(elambda, "function.name") 300 301 lmu <- as.list(substitute(lmu)) 302 emu <- link2list(lmu) 303 lmu <- attr(emu, "function.name") 304 305 lsigma <- as.list(substitute(lsigma)) 306 esigma <- link2list(lsigma) 307 lsigma <- attr(esigma, "function.name") 308 309 310 if (!is.Numeric(ilambda)) 311 stop("bad input for argument 'ilambda'") 312 if (length(isigma) && !is.Numeric(isigma, positive = TRUE)) 313 stop("bad input for argument 'isigma'") 314 315 new("vglmff", 316 blurb = c("LMS Quantile Regression ", 317 "(Box-Cox transformation to a Gamma distribution)\n", 318 "Links: ", 319 namesof("lambda", link = llambda, earg = elambda), ", ", 320 namesof("mu", link = lmu, earg = emu), ", ", 321 namesof("sigma", link = lsigma, earg = esigma)), 322 constraints = eval(substitute(expression({ 323 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 324 predictors.names = predictors.names, 325 M1 = 3) 326 }), list(.zero = zero))), 327 328 infos = eval(substitute(function(...) { 329 list(M1 = 3, 330 Q1 = 1, 331 expected = TRUE, 332 multipleResponses = FALSE, 333 parameters.names = c("lambda", "mu", "sigma"), 334 llambda = .llambda , 335 lmu = .lmu , 336 lsigma = .lsigma , 337 percentiles = .percentiles , # For the original fit only 338 true.mu = FALSE, # quantiles 339 zero = .zero ) 340 }, list( .zero = zero, 341 .percentiles = percentiles, 342 .llambda = llambda, .lmu = lmu, .lsigma = lsigma ))), 343 344 initialize = eval(substitute(expression({ 345 346 w.y.check(w = w, y = y, 347 Is.positive.y = TRUE, 348 ncol.w.max = 1, ncol.y.max = 1) 349 350 predictors.names <- c( 351 namesof("lambda", .llambda, earg = .elambda, short = TRUE), 352 namesof("mu", .lmu, earg = .emu, short = TRUE), 353 namesof("sigma", .lsigma, earg = .esigma, short = TRUE)) 354 355 extra$percentiles <- .percentiles 356 357 358 if (!length(etastart)) { 359 360 Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 361 y = y, w = w, df = .idf.mu ) 362 fv.init <- c(predict(Fit5, x = x[, min(ncol(x), 2)])$y) 363 364 lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0 365 366 sigma.init <- if (is.null( .isigma )) { 367 myratio <- ((y/fv.init)^lambda.init-1) / lambda.init 368 if (is.numeric( .idf.sigma ) && 369 is.finite( .idf.sigma )) { 370 fit600 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 371 y = (myratio)^2, 372 w = w, df = .idf.sigma ) 373 sqrt(c(abs(predict(fit600, x = x[, min(ncol(x), 2)])$y))) 374 } else { 375 sqrt(var(myratio)) 376 } 377 } else .isigma 378 379 etastart <- 380 cbind(theta2eta(lambda.init, .llambda , earg = .elambda ), 381 theta2eta(fv.init, .lmu , earg = .emu ), 382 theta2eta(sigma.init, .lsigma , earg = .esigma )) 383 } 384 }), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 385 .elambda = elambda, .emu = emu, .esigma = esigma, 386 .idf.mu = idf.mu, 387 .idf.sigma = idf.sigma, 388 .ilambda = ilambda, .isigma = isigma, 389 .percentiles = percentiles ))), 390 391 392 linkinv = eval(substitute(function(eta, extra = NULL) { 393 pcent <- extra$percentiles 394 395 eta[, 1] <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 396 eta[, 2] <- eta2theta(eta[, 2], .lmu , earg = .emu ) 397 eta[, 3] <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 398 qtplot.lms.bcg(percentiles = pcent, eta = eta) 399 }, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 400 .elambda = elambda, .emu = emu, .esigma = esigma ))), 401 last = eval(substitute(expression({ 402 misc$link <- c(lambda = .llambda , mu = .lmu , sigma = .lsigma ) 403 404 misc$earg <- list(lambda = .elambda , mu = .emu , sigma = .esigma ) 405 406 misc$percentiles <- .percentiles # These are argument values 407 if (control$cdf) { 408 post$cdf <- cdf.lms.bcg(y, eta0 = matrix(c(lambda, mymu, sigma), 409 ncol = 3, dimnames = list(dimnames(x)[[1]], NULL))) 410 } 411 }), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 412 .elambda = elambda, .emu = emu, .esigma = esigma, 413 .percentiles = percentiles ))), 414 loglikelihood = eval(substitute( 415 function(mu, y, w, residuals = FALSE, eta, 416 extra = NULL, 417 summation = TRUE) { 418 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 419 mu <- eta2theta(eta[, 2], .lmu , earg = .emu ) 420 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 421 Gee <- (y / mu)^lambda 422 theta <- 1 / (sigma * lambda)^2 423 if (residuals) { 424 stop("loglikelihood residuals not implemented yet") 425 } else { 426 ll.elts <- c(w) * (log(abs(lambda)) + theta * (log(theta) + 427 log(Gee)-Gee) - lgamma(theta) - log(y)) 428 if (summation) { 429 sum(ll.elts) 430 } else { 431 ll.elts 432 } 433 } 434 }, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 435 .elambda = elambda, .emu = emu, .esigma = esigma ))), 436 vfamily = c("lms.bcg", "lmscreg"), 437 validparams = eval(substitute(function(eta, y, extra = NULL) { 438 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 439 mymu <- eta2theta(eta[, 2], .lmu , earg = .emu ) 440 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 441 okay1 <- all(is.finite(mymu )) && 442 all(is.finite(sigma )) && all(0 < sigma) && 443 all(is.finite(lambda)) 444 okay1 445 }, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 446 .elambda = elambda, .emu = emu, .esigma = esigma ))), 447 deriv = eval(substitute(expression({ 448 lambda <- eta2theta(eta[, 1], .llambda, earg = .elambda ) 449 mymu <- eta2theta(eta[, 2], .lmu, earg = .emu ) 450 sigma <- eta2theta(eta[, 3], .lsigma, earg = .esigma ) 451 452 Gee <- (y / mymu)^lambda 453 theta <- 1 / (sigma * lambda)^2 454 dd <- digamma(theta) 455 456 dl.dlambda <- (1 + 2 * theta * (dd + Gee -1 -log(theta) - 457 0.5 * (Gee + 1) * log(Gee))) / lambda 458 dl.dmu <- lambda * theta * (Gee-1) / mymu 459 dl.dsigma <- 2*theta*(dd + Gee - log(theta * Gee)-1) / sigma 460 461 dlambda.deta <- dtheta.deta(lambda, link = .llambda , earg = .elambda ) 462 dmu.deta <- dtheta.deta(mymu, link = .lmu , earg = .emu ) 463 dsigma.deta <- dtheta.deta(sigma, link = .lsigma , earg = .esigma ) 464 465 cbind(dl.dlambda * dlambda.deta, 466 dl.dmu * dmu.deta, 467 dl.dsigma * dsigma.deta) * w 468 }), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 469 .elambda = elambda, .emu = emu, .esigma = esigma ))), 470 weight = eval(substitute(expression({ 471 tritheta <- trigamma(theta) 472 wz <- matrix(0, n, 6) 473 474 if (TRUE) { 475 part2 <- dd + 2/theta - 2*log(theta) 476 wz[,iam(1, 1, M)] <- ((1 + theta*(tritheta*(1+4*theta) - 477 4*(1+1/theta) - log(theta)*(2/theta - 478 log(theta)) + dd*part2)) / lambda^2) * 479 dlambda.deta^2 480 } else { 481 temp <- mean( Gee*(log(Gee))^2 ) 482 wz[,iam(1, 1, M)] <- ((4 * theta * (theta * tritheta-1) - 1 + 483 theta*temp) / lambda^2) * dlambda.deta^2 484 } 485 486 wz[,iam(2, 2, M)] <- dmu.deta^2 / (mymu * sigma)^2 487 wz[,iam(3, 3, M)] <- (4 * theta * (theta * tritheta - 1) / sigma^2) * 488 dsigma.deta^2 489 wz[,iam(1, 2, M)] <- (-theta * (dd + 1 / theta - log(theta)) / mymu) * 490 dlambda.deta * dmu.deta 491 wz[,iam(1, 3, M)] <- 2 * theta^1.5 * (2 * theta * tritheta - 2 - 492 1 / theta) * dlambda.deta * dsigma.deta 493 c(w) * wz 494 }), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma, 495 .elambda = elambda, .emu = emu, .esigma = esigma )))) 496} 497 498 499 500 501dy.dpsi.yeojohnson <- function(psi, lambda) { 502 503 L <- max(length(psi), length(lambda)) 504 if (length(psi) != L) psi <- rep_len(psi, L) 505 if (length(lambda) != L) lambda <- rep_len(lambda, L) 506 507 ifelse(psi > 0, (1 + psi * lambda)^(1/lambda - 1), 508 (1 - (2-lambda) * psi)^((lambda - 1) / (2-lambda))) 509} 510 511 512dyj.dy.yeojohnson <- function(y, lambda) { 513 L <- max(length(y), length(lambda)) 514 if (length(y) != L) y <- rep_len(y, L) 515 if (length(lambda) != L) lambda <- rep_len(lambda, L) 516 517 ifelse(y>0, (1 + y)^(lambda - 1), (1 - y)^(1 - lambda)) 518} 519 520 521 522 yeo.johnson <- function(y, lambda, derivative = 0, 523 epsilon = sqrt(.Machine$double.eps), 524 inverse = FALSE) { 525 526 if (!is.Numeric(derivative, length.arg = 1, 527 integer.valued = TRUE) || 528 derivative < 0) 529 stop("argument 'derivative' must be a non-negative integer") 530 531 ans <- y 532 if (!is.Numeric(epsilon, length.arg = 1, positive = TRUE)) 533 stop("argument 'epsilon' must be a single positive number") 534 L <- max(length(lambda), length(y)) 535 if (length(y) != L) y <- rep_len(y, L) 536 if (length(lambda) != L) lambda <- rep_len(lambda, L) 537 538 if (inverse) { 539 if (derivative != 0) 540 stop("argument 'derivative' must 0 when inverse = TRUE") 541 if (any(index <- y >= 0 & abs(lambda ) > epsilon)) 542 ans[index] <- (y[index]*lambda[index] + 1)^(1/lambda[index]) - 1 543 if (any(index <- y >= 0 & abs(lambda ) <= epsilon)) 544 ans[index] <- expm1(y[index]) 545 if (any(index <- y < 0 & abs(lambda-2) > epsilon)) 546 ans[index] <- 1 - (-(2-lambda[index]) * 547 y[index]+1)^(1/(2-lambda[index])) 548 if (any(index <- y < 0 & abs(lambda-2) <= epsilon)) 549 ans[index] <- -expm1(-y[index]) 550 return(ans) 551 } 552 if (derivative == 0) { 553 if (any(index <- y >= 0 & abs(lambda ) > epsilon)) 554 ans[index] <- ((y[index]+1)^(lambda[index]) - 1) / lambda[index] 555 if (any(index <- y >= 0 & abs(lambda ) <= epsilon)) 556 ans[index] <- log1p(y[index]) 557 if (any(index <- y < 0 & abs(lambda-2) > epsilon)) 558 ans[index] <- -((-y[index]+1)^(2-lambda[index]) - 1)/(2 - 559 lambda[index]) 560 if (any(index <- y < 0 & abs(lambda-2) <= epsilon)) 561 ans[index] <- -log1p(-y[index]) 562 } else { 563 psi <- Recall(y = y, lambda = lambda, derivative = derivative - 1, 564 epsilon = epsilon, inverse = inverse) 565 if (any(index <- y >= 0 & abs(lambda ) > epsilon)) 566 ans[index] <- ( (y[index]+1)^(lambda[index]) * 567 (log1p(y[index]))^(derivative) - derivative * 568 psi[index] ) / lambda[index] 569 if (any(index <- y >= 0 & abs(lambda ) <= epsilon)) 570 ans[index] <- (log1p(y[index]))^(derivative + 1) / (derivative + 1) 571 if (any(index <- y < 0 & abs(lambda-2) > epsilon)) 572 ans[index] <- -( (-y[index]+1)^(2-lambda[index]) * 573 (-log1p(-y[index]))^(derivative) - derivative * 574 psi[index] ) / (2-lambda[index]) 575 if (any(index <- y < 0 & abs(lambda-2) <= epsilon)) 576 ans[index] <- (-log1p(-y[index]))^(derivative + 1) / (derivative + 1) 577 } 578 ans 579} 580 581 582dpsi.dlambda.yjn <- function(psi, lambda, mymu, sigma, 583 derivative = 0, smallno = 1.0e-8) { 584 585 if (!is.Numeric(derivative, length.arg = 1, 586 integer.valued = TRUE) || 587 derivative < 0) 588 stop("argument 'derivative' must be a non-negative integer") 589 if (!is.Numeric(smallno, length.arg = 1, positive = TRUE)) 590 stop("argument 'smallno' must be a single positive number") 591 592 L <- max(length(psi), length(lambda), length(mymu), length(sigma)) 593 if (length(psi) != L) psi <- rep_len(psi, L) 594 if (length(lambda) != L) lambda <- rep_len(lambda, L) 595 if (length(mymu) != L) mymu <- rep_len(mymu, L) 596 if (length(sigma) != L) sigma <- rep_len(sigma, L) 597 598 answer <- matrix(NA_real_, L, derivative+1) 599 CC <- psi >= 0 600 BB <- ifelse(CC, lambda, -2+lambda) 601 AA <- psi * BB 602 temp8 <- if (derivative > 0) { 603 answer[,1:derivative] <- 604 Recall(psi = psi, lambda = lambda, mymu = mymu, sigma = sigma, 605 derivative = derivative-1, smallno = smallno) 606 answer[,derivative] * derivative 607 } else { 608 0 609 } 610 answer[, 1+derivative] <- ((AA+1) * 611 (log1p(AA)/BB)^derivative - temp8) / BB 612 613 pos <- (CC & abs(lambda) <= smallno) | (!CC & abs(lambda-2) <= smallno) 614 if (any(pos)) 615 answer[pos,1+derivative] = 616 (answer[pos, 1]^(1+derivative))/(derivative+1) 617 answer 618} 619 620 621 622gh.weight.yjn.11 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 623 624 625 if (length(derivmat)) { 626 ((derivmat[, 2]/sigma)^2 + 627 sqrt(2) * z * derivmat[, 3] / sigma) / sqrt(pi) 628 } else { 629 # Long-winded way 630 psi <- mymu + sqrt(2) * sigma * z 631 (1 / sqrt(pi)) * 632 (dpsi.dlambda.yjn(psi, lambda, mymu, sigma, 633 derivative = 1)[, 2]^2 + 634 (psi - mymu) * 635 dpsi.dlambda.yjn(psi, lambda, mymu, sigma, 636 derivative = 2)[, 3]) / sigma^2 637 } 638} 639 640 641 642gh.weight.yjn.12 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 643 if (length(derivmat)) { 644 (-derivmat[, 2]) / (sqrt(pi) * sigma^2) 645 } else { 646 psi <- mymu + sqrt(2) * sigma * z 647 (1 / sqrt(pi)) * (-dpsi.dlambda.yjn(psi, lambda, mymu, sigma, 648 derivative = 1)[, 2]) / sigma^2 649 } 650} 651 652 653 654gh.weight.yjn.13 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 655 if (length(derivmat)) { 656 sqrt(8 / pi) * (-derivmat[, 2]) * z / sigma^2 657 } else { 658 psi <- mymu + sqrt(2) * sigma * z 659 (1 / sqrt(pi)) * 660 (-2 * dpsi.dlambda.yjn(psi, lambda, mymu, sigma, 661 derivative = 1)[, 2]) * 662 (psi - mymu) / sigma^3 663 } 664} 665 666 667 668glag.weight.yjn.11 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 669 670 671 if (length(derivmat)) { 672 derivmat[, 4] * (derivmat[, 2]^2 + sqrt(2) * sigma * z * derivmat[, 3]) 673 } else { 674 psi <- mymu + sqrt(2) * sigma * z 675 discontinuity <- -mymu / (sqrt(2) * sigma) 676 (1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) * 677 (1 / sqrt(pi)) * 678 (dpsi.dlambda.yjn(psi, lambda, mymu, sigma, derivative = 1)[, 2]^2 + 679 (psi - mymu) * 680 dpsi.dlambda.yjn(psi, lambda, mymu, 681 sigma, derivative = 2)[, 3]) / sigma^2 682 } 683} 684 685 686 687glag.weight.yjn.12 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 688 discontinuity <- -mymu / (sqrt(2) * sigma) 689 if (length(derivmat)) { 690 derivmat[, 4] * (-derivmat[, 2]) 691 } else { 692 psi <- mymu + sqrt(2) * sigma * z 693 (1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) * 694 (1 / sqrt(pi)) * 695 (- dpsi.dlambda.yjn(psi, lambda, mymu, 696 sigma, derivative = 1)[, 2]) / sigma^2 697 } 698} 699 700 701 702glag.weight.yjn.13 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 703 if (length(derivmat)) { 704 derivmat[, 4] * (-derivmat[, 2]) * sqrt(8) * z 705 } else { 706 psi <- mymu + sqrt(2) * sigma * z 707 discontinuity <- -mymu / (sqrt(2) * sigma) 708 (1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) * 709 (1 / sqrt(pi)) * 710 (-2 * dpsi.dlambda.yjn(psi, lambda, mymu, 711 sigma, derivative = 1)[, 2]) * 712 (psi - mymu) / sigma^3 713 } 714} 715 716 717 718gleg.weight.yjn.11 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 719 720 721 722 723 if (length(derivmat)) { 724 derivmat[, 4] * (derivmat[, 2]^2 + sqrt(2) * sigma*z* derivmat[, 3]) 725 } else { 726 psi <- mymu + sqrt(2) * sigma * z 727 (exp(-z^2) / sqrt(pi)) * 728 (dpsi.dlambda.yjn(psi, lambda, mymu, sigma, derivative = 1)[, 2]^2 + 729 (psi - mymu) * 730 dpsi.dlambda.yjn(psi, lambda, mymu, sigma, 731 derivative = 2)[, 3]) / sigma^2 732 } 733} 734 735 736 737gleg.weight.yjn.12 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 738 if (length(derivmat)) { 739 derivmat[, 4] * (- derivmat[, 2]) 740 } else { 741 psi <- mymu + sqrt(2) * sigma * z 742 (exp(-z^2) / sqrt(pi)) * 743 (- dpsi.dlambda.yjn(psi, lambda, mymu, sigma, 744 derivative = 1)[, 2]) / sigma^2 745 } 746} 747 748 749 750gleg.weight.yjn.13 <- function(z, lambda, mymu, sigma, derivmat = NULL) { 751 if (length(derivmat)) { 752 derivmat[, 4] * (-derivmat[, 2]) * sqrt(8) * z 753 } else { 754 psi <- mymu + sqrt(2) * sigma * z 755 (exp(-z^2) / sqrt(pi)) * 756 (-2 * dpsi.dlambda.yjn(psi, lambda, mymu, 757 sigma, derivative = 1)[, 2]) * 758 (psi - mymu) / sigma^3 759 } 760} 761 762 763 764 765lms.yjn2.control <- function(save.weights = TRUE, ...) { 766 list(save.weights = save.weights) 767} 768 769 lms.yjn2 <- function(percentiles = c(25, 50, 75), 770 zero = c("lambda", "sigma"), 771 llambda = "identitylink", 772 lmu = "identitylink", 773 lsigma = "loglink", 774 idf.mu = 4, 775 idf.sigma = 2, 776 ilambda = 1.0, 777 isigma = NULL, 778 yoffset = NULL, 779 nsimEIM = 250) { 780 781 llambda <- as.list(substitute(llambda)) 782 elambda <- link2list(llambda) 783 llambda <- attr(elambda, "function.name") 784 785 lmu <- as.list(substitute(lmu)) 786 emu <- link2list(lmu) 787 lmu <- attr(emu, "function.name") 788 789 lsigma <- as.list(substitute(lsigma)) 790 esigma <- link2list(lsigma) 791 lsigma <- attr(esigma, "function.name") 792 793 794 795 if (!is.Numeric(ilambda)) 796 stop("bad input for argument 'ilambda'") 797 if (length(isigma) && 798 !is.Numeric(isigma, positive = TRUE)) 799 stop("bad input for argument 'isigma'") 800 801 new("vglmff", 802 blurb = c("LMS Quantile Regression (Yeo-Johnson transformation", 803 " to normality)\n", 804 "Links: ", 805 namesof("lambda", link = llambda, earg = elambda), ", ", 806 namesof("mu", link = lmu, earg = emu ), ", ", 807 namesof("sigma", link = lsigma, earg = esigma )), 808 constraints = eval(substitute(expression({ 809 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 810 predictors.names = predictors.names, 811 M1 = 3) 812 }), list( .zero = zero ))), 813 814 infos = eval(substitute(function(...) { 815 list(M1 = 3, 816 Q1 = 1, 817 expected = TRUE, 818 multipleResponses = FALSE, 819 parameters.names = c("lambda", "mu", "sigma"), 820 llambda = .llambda , 821 lmu = .lmu , 822 lsigma = .lsigma , 823 percentiles = .percentiles , # For the original fit only 824 true.mu = FALSE, # quantiles 825 zero = .zero ) 826 }, list( .zero = zero, 827 .percentiles = percentiles, 828 .llambda = llambda, .lmu = lmu, .lsigma = lsigma ))), 829 830 initialize = eval(substitute(expression({ 831 832 w.y.check(w = w, y = y, 833 ncol.w.max = 1, ncol.y.max = 1) 834 835 836 extra$percentiles <- .percentiles 837 838 839 predictors.names <- 840 c(namesof("lambda", .llambda, earg = .elambda, short= TRUE), 841 namesof("mu", .lmu, earg = .emu, short= TRUE), 842 namesof("sigma", .lsigma, earg = .esigma, short= TRUE)) 843 844 y.save <- y 845 yoff <- if (is.Numeric( .yoffset)) .yoffset else -median(y) 846 extra$yoffset <- yoff 847 y <- y + yoff 848 849 if (!length(etastart)) { 850 lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1. 851 852 y.tx <- yeo.johnson(y, lambda.init) 853 fv.init = 854 if (smoothok <- 855 (length(unique(sort(x[, min(ncol(x), 2)]))) > 7)) { 856 fit700 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 857 y = y.tx, w = w, df = .idf.mu ) 858 c(predict(fit700, x = x[, min(ncol(x), 2)])$y) 859 } else { 860 rep_len(weighted.mean(y, w), n) 861 } 862 863 sigma.init <- if (!is.Numeric(.isigma)) { 864 if (is.Numeric( .idf.sigma) && smoothok) { 865 fit710 = vsmooth.spline(x = x[, min(ncol(x), 2)], 866 y = (y.tx - fv.init)^2, 867 w = w, df = .idf.sigma) 868 sqrt(c(abs(predict(fit710, 869 x = x[, min(ncol(x), 2)])$y))) 870 } else { 871 sqrt( sum( w * (y.tx - fv.init)^2 ) / sum(w) ) 872 } 873 } else 874 .isigma 875 876 etastart <- matrix(0, n, 3) 877 etastart[, 1] <- theta2eta(lambda.init, .llambda, earg = .elambda) 878 etastart[, 2] <- theta2eta(fv.init, .lmu, earg = .emu) 879 etastart[, 3] <- theta2eta(sigma.init, .lsigma, earg = .esigma) 880 881 } 882 }), list(.llambda = llambda, .lmu = lmu, .lsigma = lsigma, 883 .elambda = elambda, .emu = emu, .esigma = esigma, 884 .ilambda = ilambda, .isigma = isigma, 885 .idf.mu = idf.mu, 886 .idf.sigma = idf.sigma, 887 .yoffset = yoffset, 888 .percentiles = percentiles ))), 889 890 linkinv = eval(substitute(function(eta, extra = NULL) { 891 pcent <- extra$percentiles 892 893 eta[, 1] <- eta2theta(eta[, 1], .llambda, earg = .elambda) 894 eta[, 3] <- eta2theta(eta[, 3], .lsigma, earg = .esigma) 895 qtplot.lms.yjn(percentiles = pcent, eta = eta, 896 yoffset = extra$yoff) 897 }, list( .esigma = esigma, .elambda = elambda, 898 .llambda = llambda, 899 .lsigma = lsigma ))), 900 last = eval(substitute(expression({ 901 misc$link <- c(lambda = .llambda, mu = .lmu, sigma = .lsigma) 902 misc$earg <- list(lambda = .elambda, mu = .emu, sigma = .esigma) 903 904 misc$nsimEIM <- .nsimEIM 905 misc$percentiles <- .percentiles # These are argument values 906 907 misc[["yoffset"]] <- extra$yoffset 908 909 y <- y.save # Restore back the value; to be attached to object 910 911 if (control$cdf) { 912 post$cdf <- cdf.lms.yjn(y + misc$yoffset, 913 eta0=matrix(c(lambda,mymu,sigma), 914 ncol=3, dimnames = list(dimnames(x)[[1]], NULL))) 915 } 916 }), list(.percentiles = percentiles, 917 .elambda = elambda, .emu = emu, .esigma = esigma, 918 .nsimEIM=nsimEIM, 919 .llambda = llambda, .lmu = lmu, .lsigma = lsigma ))), 920 loglikelihood = eval(substitute( 921 function(mu, y, w, residuals = FALSE, eta, 922 extra = NULL, 923 summation = TRUE) { 924 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 925 mu <- eta2theta(eta[, 2], .lmu , earg = .emu ) 926 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 927 psi <- yeo.johnson(y, lambda) 928 if (residuals) { 929 stop("loglikelihood residuals not implemented yet") 930 } else { 931 ll.elts <- c(w) * (-log(sigma) - 0.5 * ((psi-mu)/sigma)^2 + 932 (lambda-1) * sign(y) * log1p(abs(y))) 933 if (summation) { 934 sum(ll.elts) 935 } else { 936 ll.elts 937 } 938 } 939 }, list( .elambda = elambda, .emu = emu, .esigma = esigma, 940 .llambda = llambda, .lmu = lmu, 941 .lsigma = lsigma ))), 942 vfamily = c("lms.yjn2", "lmscreg"), 943 validparams = eval(substitute(function(eta, y, extra = NULL) { 944 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 945 mymu <- eta2theta(eta[, 2], .lmu , earg = .emu ) 946 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 947 okay1 <- all(is.finite(mymu )) && 948 all(is.finite(sigma )) && all(0 < sigma) && 949 all(is.finite(lambda)) 950 okay1 951 }, list( .elambda = elambda, .emu = emu, .esigma = esigma, 952 .llambda = llambda, .lmu = lmu, 953 .lsigma = lsigma ))), 954 deriv = eval(substitute(expression({ 955 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 956 mymu <- eta2theta(eta[, 2], .lmu , earg = .emu ) 957 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 958 dlambda.deta <- dtheta.deta(lambda, link = .llambda, earg = .elambda) 959 dmu.deta <- dtheta.deta(mymu, link = .lmu, earg = .emu) 960 dsigma.deta <- dtheta.deta(sigma, link = .lsigma, earg = .esigma) 961 962 psi <- yeo.johnson(y, lambda) 963 d1 <- yeo.johnson(y, lambda, deriv = 1) 964 AA <- (psi - mymu) / sigma 965 dl.dlambda <- -AA * d1 /sigma + sign(y) * log1p(abs(y)) 966 dl.dmu <- AA / sigma 967 dl.dsigma <- (AA^2 -1) / sigma 968 dthetas.detas <- cbind(dlambda.deta, dmu.deta, dsigma.deta) 969 c(w) * cbind(dl.dlambda, dl.dmu, dl.dsigma) * dthetas.detas 970 }), list( .elambda = elambda, .emu = emu, .esigma = esigma, 971 .llambda = llambda, .lmu = lmu, 972 .lsigma = lsigma ))), 973 weight = eval(substitute(expression({ 974 975 976 run.varcov <- 0 977 ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE) 978 for (ii in 1:( .nsimEIM )) { 979 psi <- rnorm(n, mymu, sigma) 980 ysim <- yeo.johnson(y = psi, lam = lambda, inverse = TRUE) 981 d1 <- yeo.johnson(ysim, lambda, deriv = 1) 982 AA <- (psi - mymu) / sigma 983 dl.dlambda <- -AA * d1 /sigma + sign(ysim) * log1p(abs(ysim)) 984 dl.dmu <- AA / sigma 985 dl.dsigma <- (AA^2 -1) / sigma 986 rm(ysim) 987 temp3 <- cbind(dl.dlambda, dl.dmu, dl.dsigma) 988 run.varcov <- ((ii-1) * run.varcov + 989 temp3[,ind1$row.index]*temp3[,ind1$col.index]) / ii 990 } 991 992 if (intercept.only) 993 run.varcov <- matrix(colMeans(run.varcov), 994 n, ncol(run.varcov), byrow = TRUE) 995 996 997 wz <- run.varcov * dthetas.detas[,ind1$row] * dthetas.detas[,ind1$col] 998 dimnames(wz) <- list(rownames(wz), NULL) # Remove the colnames 999 c(w) * wz 1000 }), list(.lsigma = lsigma, 1001 .esigma = esigma, .elambda = elambda, 1002 .nsimEIM=nsimEIM, 1003 .llambda = llambda)))) 1004} 1005 1006 1007 1008 1009 lms.yjn <- 1010 function(percentiles = c(25, 50, 75), 1011 zero = c("lambda", "sigma"), 1012 llambda = "identitylink", 1013 lsigma = "loglink", 1014 idf.mu = 4, 1015 idf.sigma = 2, 1016 ilambda = 1.0, 1017 isigma = NULL, 1018 rule = c(10, 5), 1019 yoffset = NULL, 1020 diagW = FALSE, iters.diagW = 6) { 1021 1022 1023 1024 1025 llambda <- as.list(substitute(llambda)) 1026 elambda <- link2list(llambda) 1027 llambda <- attr(elambda, "function.name") 1028 1029 lsigma <- as.list(substitute(lsigma)) 1030 esigma <- link2list(lsigma) 1031 lsigma <- attr(esigma, "function.name") 1032 1033 1034 1035 rule <- rule[1] # Nbr of points (common) for all the quadrature schemes 1036 if (rule != 5 && rule != 10) 1037 stop("only rule=5 or 10 is supported") 1038 1039 new("vglmff", 1040 blurb = c("LMS Quantile Regression ", 1041 "(Yeo-Johnson transformation to normality)\n", 1042 "Links: ", 1043 namesof("lambda", link = llambda, earg = elambda), ", ", 1044 namesof("mu", link = "identitylink", earg = list()), ", ", 1045 namesof("sigma", link = lsigma, earg = esigma)), 1046 constraints = eval(substitute(expression({ 1047 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1048 predictors.names = predictors.names, 1049 M1 = 3) 1050 }), list(.zero = zero))), 1051 1052 infos = eval(substitute(function(...) { 1053 list(M1 = 3, 1054 Q1 = 1, 1055 expected = TRUE, 1056 multipleResponses = FALSE, 1057 parameters.names = c("lambda", "mu", "sigma"), 1058 llambda = .llambda , 1059 lmu = "identitylink", 1060 lsigma = .lsigma , 1061 percentiles = .percentiles , # For the original fit only 1062 true.mu = FALSE, # quantiles 1063 zero = .zero ) 1064 }, list( .zero = zero, 1065 .percentiles = percentiles, 1066 .llambda = llambda, .lsigma = lsigma ))), 1067 1068 initialize = eval(substitute(expression({ 1069 1070 w.y.check(w = w, y = y, 1071 ncol.w.max = 1, ncol.y.max = 1) 1072 1073 1074 1075 predictors.names <- 1076 c(namesof("lambda", .llambda, earg = .elambda , short = TRUE), 1077 "mu", 1078 namesof("sigma", .lsigma, earg = .esigma , short = TRUE)) 1079 1080 extra$percentiles <- .percentiles 1081 1082 1083 y.save <- y 1084 yoff <- if (is.Numeric( .yoffset )) .yoffset else -median(y) 1085 extra$yoffset <- yoff 1086 y <- y + yoff 1087 1088 if (!length(etastart)) { 1089 1090 lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0 1091 1092 y.tx <- yeo.johnson(y, lambda.init) 1093 if (smoothok <- 1094 (length(unique(sort(x[, min(ncol(x), 2)]))) > 7)) { 1095 fit700 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 1096 y = y.tx, w = w, df = .idf.mu ) 1097 fv.init <- c(predict(fit700, x = x[, min(ncol(x), 2)])$y) 1098 } else { 1099 fv.init <- rep_len(weighted.mean(y, w), n) 1100 } 1101 1102 sigma.init <- if (!is.Numeric( .isigma )) { 1103 if (is.Numeric( .idf.sigma) && 1104 smoothok) { 1105 fit710 = vsmooth.spline(x = x[, min(ncol(x), 2)], 1106 y = (y.tx - fv.init)^2, 1107 w = w, df = .idf.sigma) 1108 sqrt(c(abs(predict(fit710, 1109 x = x[, min(ncol(x), 2)])$y))) 1110 } else { 1111 sqrt( sum( w * (y.tx - fv.init)^2 ) / sum(w) ) 1112 } 1113 } else 1114 .isigma 1115 1116 etastart <- 1117 cbind(theta2eta(lambda.init, .llambda , earg = .elambda ), 1118 fv.init, 1119 theta2eta(sigma.init, .lsigma , earg = .esigma )) 1120 1121 } 1122 }), list( .llambda = llambda, .lsigma = lsigma, 1123 .elambda = elambda, .esigma = esigma, 1124 .ilambda = ilambda, .isigma = isigma, 1125 .idf.mu = idf.mu, 1126 .idf.sigma = idf.sigma, 1127 .yoffset = yoffset, 1128 .percentiles = percentiles ))), 1129 1130 1131 linkinv = eval(substitute(function(eta, extra = NULL) { 1132 pcent <- extra$percentiles 1133 1134 eta[, 1] <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 1135 eta[, 3] <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 1136 qtplot.lms.yjn(percentiles = pcent, 1137 eta = eta, yoffset = extra$yoff) 1138 }, list(.percentiles = percentiles, 1139 .esigma = esigma, 1140 .elambda = elambda, 1141 .llambda = llambda, 1142 .lsigma = lsigma))), 1143 last = eval(substitute(expression({ 1144 misc$link <- c(lambda = .llambda , 1145 mu = "identitylink", 1146 sigma = .lsigma ) 1147 1148 misc$earg <- list(lambda = .elambda , 1149 mu = list(theta = NULL), 1150 sigma = .esigma ) 1151 1152 misc$percentiles <- .percentiles # These are argument values 1153 misc$true.mu <- FALSE # $fitted is not a true mu 1154 misc[["yoffset"]] <- extra$yoff 1155 1156 y <- y.save # Restore back the value; to be attached to object 1157 1158 if (control$cdf) { 1159 post$cdf = 1160 cdf.lms.yjn(y + misc$yoffset, 1161 eta0 = matrix(c(lambda,mymu,sigma), 1162 ncol = 3, 1163 dimnames = list(dimnames(x)[[1]], NULL))) 1164 } 1165 }), list( .percentiles = percentiles, 1166 .elambda = elambda, .esigma = esigma, 1167 .llambda = llambda, .lsigma = lsigma))), 1168 loglikelihood = eval(substitute( 1169 function(mu, y, w, residuals = FALSE, eta, 1170 extra = NULL, summation = TRUE) { 1171 1172 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 1173 mu <- eta[, 2] 1174 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 1175 psi <- yeo.johnson(y, lambda) 1176 if (residuals) { 1177 stop("loglikelihood residuals not implemented yet") 1178 } else { 1179 ll.elts <- c(w) * (-log(sigma) - 0.5 * ((psi-mu)/sigma)^2 + 1180 (lambda-1) * sign(y) * log1p(abs(y))) 1181 if (summation) { 1182 sum(ll.elts) 1183 } else { 1184 ll.elts 1185 } 1186 } 1187 }, list( .esigma = esigma, .elambda = elambda, 1188 .lsigma = lsigma, .llambda = llambda))), 1189 vfamily = c("lms.yjn", "lmscreg"), 1190 validparams = eval(substitute(function(eta, y, extra = NULL) { 1191 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 1192 mymu <- eta[, 2] 1193 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 1194 okay1 <- all(is.finite(mymu )) && 1195 all(is.finite(sigma )) && all(0 < sigma) && 1196 all(is.finite(lambda)) 1197 okay1 1198 }, list( .esigma = esigma, .elambda = elambda, 1199 .lsigma = lsigma, .llambda = llambda))), 1200 deriv = eval(substitute(expression({ 1201 lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda ) 1202 mymu <- eta[, 2] 1203 sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma ) 1204 1205 psi <- yeo.johnson(y, lambda) 1206 d1 <- yeo.johnson(y, lambda, deriv = 1) 1207 AA <- (psi - mymu) / sigma 1208 1209 dl.dlambda <- -AA * d1 /sigma + sign(y) * log1p(abs(y)) 1210 dl.dmu <- AA / sigma 1211 dl.dsigma <- (AA^2 -1) / sigma 1212 dlambda.deta <- dtheta.deta(lambda, link = .llambda, earg = .elambda ) 1213 dsigma.deta <- dtheta.deta(sigma, link = .lsigma, earg = .esigma ) 1214 1215 cbind(dl.dlambda * dlambda.deta, 1216 dl.dmu, 1217 dl.dsigma * dsigma.deta) * c(w) 1218 }), list( .esigma = esigma, .elambda = elambda, 1219 .lsigma = lsigma, .llambda = llambda ))), 1220 weight = eval(substitute(expression({ 1221 wz <- matrix(0, n, 6) 1222 1223 1224 wz[,iam(2, 2, M)] <- 1 / sigma^2 1225 wz[,iam(3, 3, M)] <- 2 * wz[,iam(2, 2, M)] # 2 / sigma^2 1226 1227 1228 if ( .rule == 10) { 1229 glag.abs = c(0.13779347054,0.729454549503, 1230 1.80834290174,3.40143369785, 1231 5.55249614006,8.33015274676, 1232 11.8437858379,16.2792578314, 1233 21.996585812, 29.9206970123) 1234 glag.wts = c(0.308441115765, 0.401119929155, 0.218068287612, 1235 0.0620874560987, 0.00950151697517, 0.000753008388588, 1236 2.82592334963e-5, 1237 4.24931398502e-7, 1.83956482398e-9, 9.91182721958e-13) 1238 } else { 1239 glag.abs = c(0.2635603197180449, 1.4134030591060496, 1240 3.5964257710396850, 1241 7.0858100058570503, 12.6408008442729685) 1242 glag.wts = c(5.217556105826727e-01, 3.986668110832433e-01, 1243 7.594244968176882e-02, 1244 3.611758679927785e-03, 2.336997238583738e-05) 1245 } 1246 1247 if ( .rule == 10) { 1248 sgh.abs = c(0.03873852801690856, 0.19823332465268367, 1249 0.46520116404433082, 1250 0.81686197962535023, 1.23454146277833154, 1251 1.70679833036403172, 1252 2.22994030591819214, 2.80910399394755972, 1253 3.46387269067033854, 1254 4.25536209637269280) 1255 sgh.wts = c(9.855210713854302e-02, 2.086780884700499e-01, 1256 2.520517066468666e-01, 1257 1.986843323208932e-01,9.719839905023238e-02, 1258 2.702440190640464e-02, 1259 3.804646170194185e-03, 2.288859354675587e-04, 1260 4.345336765471935e-06, 1261 1.247734096219375e-08) 1262 } else { 1263 sgh.abs = c(0.1002421519682381, 0.4828139660462573, 1264 1.0609498215257607, 1265 1.7797294185202606, 2.6697603560875995) 1266 sgh.wts = c(0.2484061520284881475,0.3923310666523834311, 1267 0.2114181930760276606, 1268 0.0332466603513424663, 0.0008248533445158026) 1269 } 1270 1271 if ( .rule == 10) { 1272 gleg.abs = c(-0.973906528517, -0.865063366689, -0.679409568299, 1273 -0.433395394129, -0.148874338982) 1274 gleg.abs = c(gleg.abs, rev(-gleg.abs)) 1275 gleg.wts = c(0.0666713443087, 0.149451349151, 0.219086362516, 1276 0.26926671931, 0.295524224715) 1277 gleg.wts = c(gleg.wts, rev(gleg.wts)) 1278 } else { 1279 gleg.abs = c(-0.9061798459386643,-0.5384693101056820, 0, 1280 0.5384693101056828, 0.9061798459386635) 1281 gleg.wts = c(0.2369268850561853,0.4786286704993680, 1282 0.5688888888888889, 1283 0.4786286704993661, 0.2369268850561916) 1284 } 1285 1286 1287 discontinuity = -mymu/(sqrt(2)*sigma) 1288 1289 1290 LL <- pmin(discontinuity, 0) 1291 UU <- pmax(discontinuity, 0) 1292 if (FALSE) { 1293 AA <- (UU-LL)/2 1294 for (kk in seq_along(gleg.wts)) { 1295 temp1 <- AA * gleg.wts[kk] 1296 abscissae <- (UU+LL)/2 + AA * gleg.abs[kk] 1297 psi <- mymu + sqrt(2) * sigma * abscissae 1298 temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma, 1299 derivative = 2) 1300 temp9 <- cbind(temp9, exp(-abscissae^2) / (sqrt(pi) * sigma^2)) 1301 1302 wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + temp1 * 1303 gleg.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9) 1304 wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + temp1 * 1305 gleg.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9) 1306 wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + temp1 * 1307 gleg.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9) 1308 } 1309 } else { 1310 temp9 <- .Fortran("yjngintf", as.double(LL), 1311 as.double(UU), 1312 as.double(gleg.abs), as.double(gleg.wts), as.integer(n), 1313 as.integer(length(gleg.abs)), as.double(lambda), 1314 as.double(mymu), as.double(sigma), answer = double(3*n), 1315 eps = as.double(1.0e-5))$ans 1316 dim(temp9) <- c(3, n) 1317 wz[,iam(1, 1, M)] <- temp9[1,] 1318 wz[,iam(1, 2, M)] <- temp9[2,] 1319 wz[,iam(1, 3, M)] <- temp9[3,] 1320 } 1321 1322 1323 1324 for (kk in seq_along(sgh.wts)) { 1325 1326 abscissae <- sign(-discontinuity) * sgh.abs[kk] 1327 psi <- mymu + sqrt(2) * sigma * abscissae # abscissae = z 1328 temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma, 1329 derivative = 2) 1330 wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + sgh.wts[kk] * 1331 gh.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9) 1332 wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + sgh.wts[kk] * 1333 gh.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9) 1334 wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + sgh.wts[kk] * 1335 gh.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9) 1336 } 1337 1338 temp1 <- exp(-discontinuity^2) 1339 for (kk in seq_along(glag.wts)) { 1340 abscissae <- sign(discontinuity) * sqrt(glag.abs[kk]) + 1341 discontinuity^2 1342 psi <- mymu + sqrt(2) * sigma * abscissae 1343 temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma, derivative = 2) 1344 temp9 <- cbind(temp9, 1345 1 / (2 * sqrt((abscissae-discontinuity^2)^2 + 1346 discontinuity^2) * 1347 sqrt(pi) * sigma^2)) 1348 temp7 <- temp1 * glag.wts[kk] 1349 wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + temp7 * 1350 glag.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9) 1351 wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + temp7 * 1352 glag.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9) 1353 wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + temp7 * 1354 glag.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9) 1355 } 1356 1357 wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] * dlambda.deta^2 1358 wz[, iam(1, 2, M)] <- wz[, iam(1, 2, M)] * dlambda.deta 1359 wz[, iam(1, 3, M)] <- wz[, iam(1, 3, M)] * dsigma.deta * dlambda.deta 1360 if ( .diagW && iter <= .iters.diagW) { 1361 wz[,iam(1, 2, M)] <- wz[, iam(1, 3, M)] <- 0 1362 } 1363 wz[, iam(2, 3, M)] <- wz[, iam(2, 3, M)] * dsigma.deta 1364 wz[, iam(3, 3, M)] <- wz[, iam(3, 3, M)] * dsigma.deta^2 1365 c(w) * wz 1366 }), list( .lsigma = lsigma, .llambda = llambda, 1367 .esigma = esigma, .elambda = elambda, 1368 .rule = rule, .diagW = diagW, 1369 .iters.diagW = iters.diagW )))) 1370} # lms.yjn 1371 1372 1373 1374 1375lmscreg.control <- function(cdf = TRUE, at.arg = NULL, x0 = NULL, ...) { 1376 1377 if (!is.logical(cdf)) { 1378 warning("'cdf' is not logical; using TRUE instead") 1379 cdf <- TRUE 1380 } 1381 list(cdf = cdf, at.arg = at.arg, x0 = x0) 1382} 1383 1384 1385 1386 1387 1388Wr1 <- function(r, w) ifelse(r <= 0, 1, w) 1389 1390 1391Wr2 <- function(r, w) (r <= 0) * 1 + (r > 0) * w 1392 1393 1394amlnormal.deviance <- function(mu, y, w, residuals = FALSE, 1395 eta, extra = NULL) { 1396 1397 M <- length(extra$w.aml) 1398 1399 if (M > 1) y <- matrix(y, extra$n, extra$M) 1400 1401 devi <- cbind((y - mu)^2) 1402 if (residuals) { 1403 stop("not sure here") 1404 wz <- VGAM.weights.function(w = w, M = extra$M, n = extra$n) 1405 return((y - mu) * sqrt(wz) * matrix(extra$w.aml,extra$n,extra$M)) 1406 } else { 1407 all.deviances <- numeric(M) 1408 myresid <- matrix(y,extra$n,extra$M) - cbind(mu) 1409 for (ii in 1:M) 1410 all.deviances[ii] <- sum(c(w) * devi[, ii] * 1411 Wr1(myresid[, ii], 1412 w = extra$w.aml[ii])) 1413 } 1414 if (is.logical(extra$individual) && extra$individual) 1415 all.deviances else sum(all.deviances) 1416} 1417 1418 1419 1420 amlnormal <- function(w.aml = 1, parallel = FALSE, 1421 lexpectile = "identitylink", 1422 iexpectile = NULL, 1423 imethod = 1, digw = 4) { 1424 1425 1426 1427 1428 if (!is.Numeric(w.aml, positive = TRUE)) 1429 stop("argument 'w.aml' must be a vector of positive values") 1430 if (!is.Numeric(imethod, length.arg = 1, 1431 integer.valued = TRUE, positive = TRUE) || 1432 imethod > 3) 1433 stop("argument 'imethod' must be 1, 2 or 3") 1434 1435 1436 1437 lexpectile <- as.list(substitute(lexpectile)) 1438 eexpectile <- link2list(lexpectile) 1439 lexpectile <- attr(eexpectile, "function.name") 1440 1441 1442 if (length(iexpectile) && !is.Numeric(iexpectile)) 1443 stop("bad input for argument 'iexpectile'") 1444 1445 new("vglmff", 1446 blurb = c("Asymmetric least squares quantile regression\n\n", 1447 "Links: ", 1448 namesof("expectile", link = lexpectile, earg = eexpectile)), 1449 constraints = eval(substitute(expression({ 1450 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1451 bool = .parallel , 1452 constraints = constraints) 1453 }), list( .parallel = parallel ))), 1454 deviance = function(mu, y, w, residuals = FALSE, eta, 1455 extra = NULL) { 1456 amlnormal.deviance(mu = mu, y = y, w = w, residuals = residuals, 1457 eta = eta, extra = extra) 1458 }, 1459 initialize = eval(substitute(expression({ 1460 extra$w.aml <- .w.aml 1461 1462 temp5 <- 1463 w.y.check(w = w, y = y, 1464 ncol.w.max = 1, ncol.y.max = 1, 1465 out.wy = TRUE, 1466 maximize = TRUE) 1467 w <- temp5$w 1468 y <- temp5$y 1469 1470 1471 1472 1473 extra$M <- M <- length(extra$w.aml) # Recycle if necessary 1474 extra$n <- n 1475 extra$y.names <- y.names <- 1476 paste("w.aml = ", round(extra$w.aml, digits = .digw ), sep = "") 1477 1478 predictors.names <- c(namesof( 1479 paste("expectile(",y.names,")", sep = ""), .lexpectile , 1480 earg = .eexpectile, tag = FALSE)) 1481 1482 if (!length(etastart)) { 1483 mean.init <- 1484 if ( .imethod == 1) 1485 rep_len(median(y), n) else 1486 if ( .imethod == 2 || .imethod == 3) 1487 rep_len(weighted.mean(y, w), n) else { 1488 junk <- lm.wfit(x = x, y = c(y), w = c(w)) 1489 junk$fitted 1490 } 1491 1492 1493 if ( .imethod == 3) 1494 mean.init <- abs(mean.init) + 0.01 1495 1496 1497 if (length( .iexpectile)) 1498 mean.init <- matrix( .iexpectile, n, M, byrow = TRUE) 1499 etastart <- 1500 matrix(theta2eta(mean.init, .lexpectile, 1501 earg = .eexpectile), n, M) 1502 } 1503 }), list( .lexpectile = lexpectile, .eexpectile = eexpectile, 1504 .iexpectile = iexpectile, 1505 .imethod = imethod, .digw = digw, .w.aml = w.aml ))), 1506 linkinv = eval(substitute(function(eta, extra = NULL) { 1507 ans <- eta <- as.matrix(eta) 1508 for (ii in 1:ncol(eta)) 1509 ans[, ii] <- eta2theta(eta[, ii], .lexpectile , earg = .eexpectile ) 1510 dimnames(ans) <- list(dimnames(eta)[[1]], extra$y.names) 1511 ans 1512 }, list( .lexpectile = lexpectile, .eexpectile = eexpectile ))), 1513 last = eval(substitute(expression({ 1514 misc$link <- rep_len(.lexpectile , M) 1515 names(misc$link) <- extra$y.names 1516 1517 misc$earg <- vector("list", M) 1518 for (ilocal in 1:M) 1519 misc$earg[[ilocal]] <- list(theta = NULL) 1520 names(misc$earg) <- names(misc$link) 1521 1522 misc$parallel <- .parallel 1523 misc$expected <- TRUE 1524 extra$percentile <- numeric(M) # These are estimates (empirical) 1525 misc$multipleResponses <- TRUE 1526 1527 1528 for (ii in 1:M) { 1529 use.w <- if (M > 1 && NCOL(w) == M) w[, ii] else w 1530 extra$percentile[ii] <- 100 * 1531 weighted.mean(myresid[, ii] <= 0, use.w) 1532 } 1533 names(extra$percentile) <- names(misc$link) 1534 1535 extra$individual <- TRUE 1536 if (!(M > 1 && NCOL(w) == M)) { 1537 extra$deviance <- 1538 amlnormal.deviance(mu = mu, y = y, w = w, 1539 residuals = FALSE, eta = eta, extra = extra) 1540 names(extra$deviance) <- extra$y.names 1541 } 1542 }), list( .lexpectile = lexpectile, 1543 .eexpectile = eexpectile, .parallel = parallel ))), 1544 vfamily = c("amlnormal"), 1545 validparams = eval(substitute(function(eta, y, extra = NULL) { 1546 mymu <- eta2theta(eta, .lexpectile , earg = .eexpectile ) 1547 okay1 <- all(is.finite(mymu)) 1548 okay1 1549 }, list( .lexpectile = lexpectile, 1550 .eexpectile = eexpectile ))), 1551 1552 deriv = eval(substitute(expression({ 1553 mymu <- eta2theta(eta, .lexpectile , earg = .eexpectile ) 1554 dexpectile.deta <- dtheta.deta(mymu, .lexpectile , 1555 earg = .eexpectile ) 1556 myresid <- matrix(y,extra$n,extra$M) - cbind(mu) 1557 wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M, 1558 byrow = TRUE)) 1559 c(w) * myresid * wor1 * dexpectile.deta 1560 }), list( .lexpectile = lexpectile, 1561 .eexpectile = eexpectile ))), 1562 1563 weight = eval(substitute(expression({ 1564 wz <- c(w) * wor1 * dexpectile.deta^2 1565 wz 1566 }), list( .lexpectile = lexpectile, 1567 .eexpectile = eexpectile )))) 1568} 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579amlpoisson.deviance <- function(mu, y, w, residuals = FALSE, eta, 1580 extra = NULL) { 1581 1582 M <- length(extra$w.aml) 1583 1584 if (M > 1) y <- matrix(y,extra$n,extra$M) 1585 1586 nz <- y > 0 1587 devi <- cbind(-(y - mu)) 1588 devi[nz] <- devi[nz] + y[nz] * log(y[nz]/mu[nz]) 1589 if (residuals) { 1590 stop("not sure here") 1591 return(sign(y - mu) * sqrt(2 * abs(devi) * w) * 1592 matrix(extra$w,extra$n,extra$M)) 1593 } else { 1594 all.deviances <- numeric(M) 1595 myresid <- matrix(y,extra$n,extra$M) - cbind(mu) 1596 for (ii in 1:M) all.deviances[ii] <- 2 * sum(c(w) * devi[, ii] * 1597 Wr1(myresid[, ii], w=extra$w.aml[ii])) 1598 } 1599 if (is.logical(extra$individual) && extra$individual) 1600 all.deviances else sum(all.deviances) 1601} 1602 1603 1604 1605 amlpoisson <- function(w.aml = 1, parallel = FALSE, imethod = 1, 1606 digw = 4, link = "loglink") { 1607 if (!is.Numeric(w.aml, positive = TRUE)) 1608 stop("'w.aml' must be a vector of positive values") 1609 1610 1611 link <- as.list(substitute(link)) 1612 earg <- link2list(link) 1613 link <- attr(earg, "function.name") 1614 1615 1616 new("vglmff", 1617 blurb = c("Poisson expectile regression by", 1618 " asymmetric maximum likelihood estimation\n\n", 1619 "Link: ", namesof("expectile", link, earg = earg)), 1620 constraints = eval(substitute(expression({ 1621 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1622 bool = .parallel , 1623 constraints = constraints) 1624 }), list( .parallel = parallel ))), 1625 deviance = function(mu, y, w, residuals = FALSE, eta, extra = NULL) { 1626 amlpoisson.deviance(mu = mu, y = y, w = w, residuals = residuals, 1627 eta = eta, extra = extra) 1628 }, 1629 initialize = eval(substitute(expression({ 1630 extra$w.aml <- .w.aml 1631 1632 temp5 <- 1633 w.y.check(w = w, y = y, 1634 ncol.w.max = 1, ncol.y.max = 1, 1635 out.wy = TRUE, 1636 maximize = TRUE) 1637 w <- temp5$w 1638 y <- temp5$y 1639 1640 1641 1642 extra$M <- M <- length(extra$w.aml) # Recycle if necessary 1643 extra$n <- n 1644 extra$y.names <- y.names <- 1645 paste("w.aml = ", round(extra$w.aml, digits = .digw ), sep = "") 1646 extra$individual <- FALSE 1647 predictors.names <- 1648 c(namesof(paste("expectile(",y.names,")", sep = ""), 1649 .link , earg = .earg , tag = FALSE)) 1650 1651 if (!length(etastart)) { 1652 mean.init <- if ( .imethod == 2) 1653 rep_len(median(y), n) else 1654 if ( .imethod == 1) 1655 rep_len(weighted.mean(y, w), n) else { 1656 junk = lm.wfit(x = x, y = c(y), w = c(w)) 1657 abs(junk$fitted) 1658 } 1659 etastart <- 1660 matrix(theta2eta(mean.init, .link , earg = .earg ), n, M) 1661 } 1662 }), list( .link = link, .earg = earg, .imethod = imethod, 1663 .digw = digw, .w.aml = w.aml ))), 1664 linkinv = eval(substitute(function(eta, extra = NULL) { 1665 mu.ans <- eta <- as.matrix(eta) 1666 for (ii in 1:ncol(eta)) 1667 mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg ) 1668 dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names) 1669 mu.ans 1670 }, list( .link = link, .earg = earg ))), 1671 last = eval(substitute(expression({ 1672 misc$multipleResponses <- TRUE 1673 misc$expected <- TRUE 1674 misc$parallel <- .parallel 1675 1676 1677 misc$link <- rep_len( .link , M) 1678 names(misc$link) <- extra$y.names 1679 1680 misc$earg <- vector("list", M) 1681 for (ilocal in 1:M) 1682 misc$earg[[ilocal]] <- list(theta = NULL) 1683 names(misc$earg) <- names(misc$link) 1684 1685 extra$percentile <- numeric(M) # These are estimates (empirical) 1686 for (ii in 1:M) 1687 extra$percentile[ii] <- 100 * weighted.mean(myresid[, ii] <= 0, w) 1688 names(extra$percentile) <- names(misc$link) 1689 1690 extra$individual <- TRUE 1691 extra$deviance <- amlpoisson.deviance(mu = mu, y = y, w = w, 1692 residuals = FALSE, 1693 eta = eta, extra = extra) 1694 names(extra$deviance) <- extra$y.names 1695 }), list( .link = link, .earg = earg, .parallel = parallel ))), 1696 linkfun = eval(substitute(function(mu, extra = NULL) { 1697 theta2eta(mu, link = .link , earg = .earg ) 1698 }, list( .link = link, .earg = earg ))), 1699 vfamily = c("amlpoisson"), 1700 validparams = eval(substitute(function(eta, y, extra = NULL) { 1701 mymu <- eta2theta(eta, .link , earg = .earg ) 1702 okay1 <- all(is.finite(mymu)) && all(0 < mymu) 1703 okay1 1704 }, list( .link = link, .earg = earg ))), 1705 deriv = eval(substitute(expression({ 1706 mymu <- eta2theta(eta, .link , earg = .earg ) 1707 dexpectile.deta <- dtheta.deta(mymu, .link , earg = .earg ) 1708 myresid <- matrix(y,extra$n,extra$M) - cbind(mu) 1709 wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M, 1710 byrow = TRUE)) 1711 c(w) * myresid * wor1 * (dexpectile.deta / mymu) 1712 }), list( .link = link, .earg = earg ))), 1713 weight = eval(substitute(expression({ 1714 use.mu <- mymu 1715 use.mu[use.mu < .Machine$double.eps^(3/4)] <- 1716 .Machine$double.eps^(3/4) 1717 wz <- c(w) * wor1 * use.mu * (dexpectile.deta / mymu)^2 1718 wz 1719 }), list( .link = link, .earg = earg )))) 1720} 1721 1722 1723 1724 1725 1726amlbinomial.deviance <- function(mu, y, w, residuals = FALSE, 1727 eta, extra = NULL) { 1728 1729 M <- length(extra$w.aml) 1730 1731 1732 if (M > 1) y <- matrix(y,extra$n,extra$M) 1733 1734 1735 devy <- y 1736 nz <- y != 0 1737 devy[nz] <- y[nz] * log(y[nz]) 1738 nz <- (1 - y) != 0 1739 devy[nz] <- devy[nz] + (1 - y[nz]) * log1p(-y[nz]) 1740 devmu <- y * log(mu) + (1 - y) * log1p(-mu) 1741 if (any(small <- mu * (1 - mu) < .Machine$double.eps)) { 1742 warning("fitted values close to 0 or 1") 1743 smu <- mu[small] 1744 sy <- y[small] 1745 smu <- ifelse(smu < .Machine$double.eps, 1746 .Machine$double.eps, smu) 1747 onemsmu <- ifelse((1 - smu) < .Machine$double.eps, 1748 .Machine$double.eps, 1 - smu) 1749 devmu[small] <- sy * log(smu) + (1 - sy) * log(onemsmu) 1750 } 1751 devi <- 2 * (devy - devmu) 1752 if (residuals) { 1753 stop("not sure here") 1754 return(sign(y - mu) * sqrt(abs(devi) * w)) 1755 } else { 1756 all.deviances <- numeric(M) 1757 myresid <- matrix(y,extra$n,extra$M) - matrix(mu,extra$n,extra$M) 1758 for (ii in 1:M) all.deviances[ii] <- sum(c(w) * devi[, ii] * 1759 Wr1(myresid[, ii], w=extra$w.aml[ii])) 1760 } 1761 if (is.logical(extra$individual) && extra$individual) 1762 all.deviances else sum(all.deviances) 1763} 1764 1765 1766 amlbinomial <- function(w.aml = 1, parallel = FALSE, digw = 4, 1767 link = "logitlink") { 1768 1769 if (!is.Numeric(w.aml, positive = TRUE)) 1770 stop("'w.aml' must be a vector of positive values") 1771 1772 1773 link <- as.list(substitute(link)) 1774 earg <- link2list(link) 1775 link <- attr(earg, "function.name") 1776 1777 1778 new("vglmff", 1779 blurb = c("Logistic expectile regression by ", 1780 "asymmetric maximum likelihood estimation\n\n", 1781 "Link: ", namesof("expectile", link, earg = earg)), 1782 constraints = eval(substitute(expression({ 1783 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1784 bool = .parallel , 1785 constraints = constraints) 1786 }), list( .parallel = parallel ))), 1787 deviance = function(mu, y, w, residuals = FALSE, eta, extra = NULL) { 1788 amlbinomial.deviance(mu = mu, y = y, w = w, residuals = residuals, 1789 eta = eta, extra = extra) 1790 }, 1791 initialize = eval(substitute(expression({ 1792 1793 1794 { 1795 NCOL <- function (x) 1796 if (is.array(x) && length(dim(x)) > 1 || 1797 is.data.frame(x)) ncol(x) else as.integer(1) 1798 1799 if (NCOL(y) == 1) { 1800 if (is.factor(y)) y <- y != levels(y)[1] 1801 nn <- rep_len(1, n) 1802 if (!all(y >= 0 & y <= 1)) 1803 stop("response values must be in [0, 1]") 1804 if (!length(mustart) && !length(etastart)) 1805 mustart <- (0.5 + w * y) / (1 + w) 1806 no.successes <- w * y 1807 if (any(abs(no.successes - round(no.successes)) > 0.001)) 1808 stop("Number of successes must be integer-valued") 1809 } else if (NCOL(y) == 2) { 1810 if (any(abs(y - round(y)) > 0.001)) 1811 stop("Count data must be integer-valued") 1812 nn <- y[, 1] + y[, 2] 1813 y <- ifelse(nn > 0, y[, 1]/nn, 0) 1814 w <- w * nn 1815 if (!length(mustart) && !length(etastart)) 1816 mustart <- (0.5 + nn * y) / (1 + nn) 1817 } else 1818 stop("Response not of the right form") 1819 } 1820 1821 extra$w.aml <- .w.aml 1822 if (ncol(y <- cbind(y)) != 1) 1823 stop("response must be a vector or a one-column matrix") 1824 extra$M <- M <- length(extra$w.aml) # Recycle if necessary 1825 extra$n <- n 1826 extra$y.names <- y.names <- 1827 paste("w.aml = ", round(extra$w.aml, digits = .digw ), 1828 sep = "") 1829 extra$individual <- FALSE 1830 predictors.names <- 1831 c(namesof(paste("expectile(", y.names, ")", sep = ""), 1832 .link , earg = .earg , tag = FALSE)) 1833 1834 if (!length(etastart)) { 1835 etastart <- matrix(theta2eta(mustart, .link , earg = .earg ), 1836 n, M) 1837 mustart <- NULL 1838 } 1839 1840 1841 }), list( .link = link, .earg = earg, 1842 .digw = digw, .w.aml = w.aml ))), 1843 linkinv = eval(substitute(function(eta, extra = NULL) { 1844 mu.ans <- eta <- as.matrix(eta) 1845 for (ii in 1:ncol(eta)) 1846 mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg ) 1847 dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names) 1848 mu.ans 1849 }, list( .link = link, .earg = earg ))), 1850 last = eval(substitute(expression({ 1851 misc$link <- rep_len(.link , M) 1852 names(misc$link) <- extra$y.names 1853 1854 misc$earg <- vector("list", M) 1855 for (ilocal in 1:M) 1856 misc$earg[[ilocal]] <- list(theta = NULL) 1857 names(misc$earg) <- names(misc$link) 1858 1859 misc$parallel <- .parallel 1860 misc$expected <- TRUE 1861 1862 extra$percentile <- numeric(M) # These are estimates (empirical) 1863 for (ii in 1:M) 1864 extra$percentile[ii] <- 100 * weighted.mean(myresid[, ii] <= 0, w) 1865 names(extra$percentile) <- names(misc$link) 1866 1867 extra$individual <- TRUE 1868 extra$deviance <- amlbinomial.deviance(mu = mu, y = y, w = w, 1869 residuals = FALSE, eta = eta, extra = extra) 1870 names(extra$deviance) <- extra$y.names 1871 }), list( .link = link, .earg = earg, .parallel = parallel ))), 1872 linkfun = eval(substitute(function(mu, extra = NULL) { 1873 theta2eta(mu, link = .link , earg = .earg ) 1874 }, list( .link = link, .earg = earg ))), 1875 vfamily = c("amlbinomial"), 1876 validparams = eval(substitute(function(eta, y, extra = NULL) { 1877 mymu <- eta2theta(eta, .link , earg = .earg ) 1878 okay1 <- all(is.finite(mymu)) && all(0 < mymu & mymu < 1) 1879 okay1 1880 }, list( .link = link, .earg = earg ))), 1881 deriv = eval(substitute(expression({ 1882 mymu <- eta2theta(eta, .link , earg = .earg ) 1883 use.mu <- mymu 1884 use.mu[use.mu < .Machine$double.eps^(3/4)] = .Machine$double.eps^(3/4) 1885 dexpectile.deta <- dtheta.deta(use.mu, .link , earg = .earg ) 1886 myresid <- matrix(y,extra$n,extra$M) - cbind(mu) 1887 wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M, 1888 byrow = TRUE)) 1889 c(w) * myresid * wor1 * (dexpectile.deta / (use.mu * (1-use.mu))) 1890 }), list( .link = link, .earg = earg ))), 1891 weight = eval(substitute(expression({ 1892 wz <- c(w) * wor1 * (dexpectile.deta^2 / (use.mu * (1 - use.mu))) 1893 wz 1894 }), list( .link = link, .earg = earg)))) 1895} 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906amlexponential.deviance <- function(mu, y, w, residuals = FALSE, 1907 eta, extra = NULL) { 1908 1909 M <- length(extra$w.aml) 1910 1911 if (M > 1) y <- matrix(y,extra$n,extra$M) 1912 1913 devy <- cbind(-log(y) - 1) 1914 devi <- cbind(-log(mu) - y / mu) 1915 if (residuals) { 1916 stop("not sure here") 1917 return(sign(y - mu) * sqrt(2 * abs(devi) * w) * 1918 matrix(extra$w,extra$n,extra$M)) 1919 } else { 1920 all.deviances <- numeric(M) 1921 myresid <- matrix(y,extra$n,extra$M) - cbind(mu) 1922 for (ii in 1:M) all.deviances[ii] = 2 * sum(c(w) * 1923 (devy[, ii] - devi[, ii]) * 1924 Wr1(myresid[, ii], w=extra$w.aml[ii])) 1925 } 1926 if (is.logical(extra$individual) && extra$individual) 1927 all.deviances else sum(all.deviances) 1928} 1929 1930 1931 1932 1933 amlexponential <- function(w.aml = 1, parallel = FALSE, imethod = 1, 1934 digw = 4, link = "loglink") { 1935 if (!is.Numeric(w.aml, positive = TRUE)) 1936 stop("'w.aml' must be a vector of positive values") 1937 if (!is.Numeric(imethod, length.arg = 1, 1938 integer.valued = TRUE, positive = TRUE) || 1939 imethod > 3) 1940 stop("argument 'imethod' must be 1, 2 or 3") 1941 1942 1943 link <- as.list(substitute(link)) 1944 earg <- link2list(link) 1945 link <- attr(earg, "function.name") 1946 1947 1948 y.names <- paste("w.aml = ", round(w.aml, digits = digw), sep = "") 1949 predictors.names <- c(namesof( 1950 paste("expectile(", y.names,")", sep = ""), link, earg = earg)) 1951 predictors.names <- paste(predictors.names, collapse = ", ") 1952 1953 1954 new("vglmff", 1955 blurb = c("Exponential expectile regression by", 1956 " asymmetric maximum likelihood estimation\n\n", 1957 "Link: ", predictors.names), 1958 constraints = eval(substitute(expression({ 1959 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1960 bool = .parallel , 1961 constraints = constraints) 1962 }), list( .parallel = parallel ))), 1963 deviance = function(mu, y, w, residuals = FALSE, eta, extra = NULL) { 1964 amlexponential.deviance(mu = mu, y = y, w = w, 1965 residuals = residuals, 1966 eta = eta, extra = extra) 1967 }, 1968 initialize = eval(substitute(expression({ 1969 extra$w.aml <- .w.aml 1970 1971 1972 temp5 <- 1973 w.y.check(w = w, y = y, 1974 Is.positive.y = TRUE, 1975 ncol.w.max = 1, ncol.y.max = 1, 1976 out.wy = TRUE, 1977 maximize = TRUE) 1978 w <- temp5$w 1979 y <- temp5$y 1980 1981 1982 1983 extra$M <- M <- length(extra$w.aml) # Recycle if necessary 1984 extra$n <- n 1985 extra$y.names <- y.names <- 1986 paste("w.aml = ", round(extra$w.aml, digits = .digw ), sep = "") 1987 extra$individual = FALSE 1988 1989 1990 predictors.names <- c(namesof( 1991 paste("expectile(", y.names, ")", sep = ""), 1992 .link , earg = .earg , tag = FALSE)) 1993 1994 if (!length(etastart)) { 1995 mean.init <- if ( .imethod == 1) 1996 rep_len(median(y), n) else 1997 if ( .imethod == 2) 1998 rep_len(weighted.mean(y, w), n) else { 1999 1 / (y + 1) 2000 } 2001 etastart <- matrix(theta2eta(mean.init, .link , earg = .earg ), 2002 n, M) 2003 } 2004 }), list( .link = link, .earg = earg, .imethod = imethod, 2005 .digw = digw, .w.aml = w.aml ))), 2006 linkinv = eval(substitute(function(eta, extra = NULL) { 2007 mu.ans <- eta <- as.matrix(eta) 2008 for (ii in 1:ncol(eta)) 2009 mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg ) 2010 dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names) 2011 mu.ans 2012 }, list( .link = link, .earg = earg ))), 2013 last = eval(substitute(expression({ 2014 misc$multipleResponses <- TRUE 2015 misc$expected <- TRUE 2016 misc$parallel <- .parallel 2017 2018 misc$link <- rep_len( .link , M) 2019 names(misc$link) <- extra$y.names 2020 2021 misc$earg <- vector("list", M) 2022 for (ilocal in 1:M) 2023 misc$earg[[ilocal]] <- list(theta = NULL) 2024 names(misc$earg) <- names(misc$link) 2025 2026 2027 extra$percentile <- numeric(M) # These are estimates (empirical) 2028 for (ii in 1:M) 2029 extra$percentile[ii] <- 100 * weighted.mean(myresid[, ii] <= 0, w) 2030 names(extra$percentile) <- names(misc$link) 2031 2032 extra$individual <- TRUE 2033 extra$deviance = 2034 amlexponential.deviance(mu = mu, y = y, w = w, 2035 residuals = FALSE, eta = eta, extra = extra) 2036 names(extra$deviance) <- extra$y.names 2037 }), list( .link = link, .earg = earg, .parallel = parallel ))), 2038 linkfun = eval(substitute(function(mu, extra = NULL) { 2039 theta2eta(mu, link = .link , earg = .earg ) 2040 }, list( .link = link, .earg = earg ))), 2041 vfamily = c("amlexponential"), 2042 validparams = eval(substitute(function(eta, y, extra = NULL) { 2043 mymu <- eta2theta(eta, .link , earg = .earg ) 2044 okay1 <- all(is.finite(mymu)) && all(0 < mymu) 2045 okay1 2046 }, list( .link = link, .earg = earg ))), 2047 deriv = eval(substitute(expression({ 2048 mymu <- eta2theta(eta, .link , earg = .earg ) 2049 bigy <- matrix(y,extra$n,extra$M) 2050 dl.dmu <- (bigy - mymu) / mymu^2 2051 2052 dmu.deta <- dtheta.deta(mymu, .link , earg = .earg ) 2053 myresid <- bigy - cbind(mymu) 2054 wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M, 2055 byrow = TRUE)) 2056 c(w) * wor1 * dl.dmu * dmu.deta 2057 }), list( .link = link, .earg = earg ))), 2058 weight = eval(substitute(expression({ 2059 ned2l.dmu2 <- 1 / mymu^2 2060 wz <- c(w) * wor1 * ned2l.dmu2 * dmu.deta^2 2061 wz 2062 }), list( .link = link, .earg = earg )))) 2063} 2064 2065 2066 2067 2068 2069 2070rho1check <- function(u, tau = 0.5) 2071 u * (tau - (u <= 0)) 2072 2073 2074 2075 2076dalap <- function(x, location = 0, scale = 1, tau = 0.5, 2077 kappa = sqrt(tau/(1-tau)), log = FALSE) { 2078 if (!is.logical(log.arg <- log) || length(log) != 1) 2079 stop("bad input for argument 'log'") 2080 rm(log) 2081 2082 2083 2084 2085 NN <- max(length(x), length(location), length(scale), length(kappa), 2086 length(tau)) 2087 if (length(x) != NN) x <- rep_len(x, NN) 2088 if (length(location) != NN) location <- rep_len(location, NN) 2089 if (length(scale) != NN) scale <- rep_len(scale, NN) 2090 if (length(kappa) != NN) kappa <- rep_len(kappa, NN) 2091 if (length(tau) != NN) tau <- rep_len(tau, NN) 2092 2093 logconst <- 0.5 * log(2) - log(scale) + log(kappa) - log1p(kappa^2) 2094 exponent <- -(sqrt(2) / scale) * abs(x - location) * 2095 ifelse(x >= location, kappa, 1/kappa) 2096 2097 indexTF <- (scale > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2098 logconst[!indexTF] <- NaN 2099 2100 if (log.arg) logconst + exponent else exp(logconst + exponent) 2101} 2102 2103 2104ralap <- function(n, location = 0, scale = 1, tau = 0.5, 2105 kappa = sqrt(tau/(1-tau))) { 2106 use.n <- if ((length.n <- length(n)) > 1) length.n else 2107 if (!is.Numeric(n, integer.valued = TRUE, 2108 length.arg = 1, positive = TRUE)) 2109 stop("bad input for argument 'n'") else n 2110 2111 location <- rep_len(location, use.n) 2112 scale <- rep_len(scale, use.n) 2113 tau <- rep_len(tau, use.n) 2114 kappa <- rep_len(kappa, use.n) 2115 ans <- location + scale * 2116 log(runif(use.n)^kappa / runif(use.n)^(1/kappa)) / sqrt(2) 2117 indexTF <- (scale > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2118 ans[!indexTF] <- NaN 2119 ans 2120} 2121 2122 2123 2124palap <- function(q, location = 0, scale = 1, tau = 0.5, 2125 kappa = sqrt(tau/(1-tau)), 2126 lower.tail = TRUE, log.p = FALSE) { 2127 2128 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 2129 stop("bad input for argument 'lower.tail'") 2130 2131 if (!is.logical(log.p) || length(log.p) != 1) 2132 stop("bad input for argument 'log.p'") 2133 2134 NN <- max(length(q), length(location), length(scale), length(kappa), 2135 length(tau)) 2136 if (length(q) != NN) q <- rep_len(q, NN) 2137 if (length(location) != NN) location <- rep_len(location, NN) 2138 if (length(scale) != NN) scale <- rep_len(scale, NN) 2139 if (length(kappa) != NN) kappa <- rep_len(kappa, NN) 2140 if (length(tau) != NN) tau <- rep_len(tau, NN) 2141 2142 exponent <- -(sqrt(2) / scale) * abs(q - location) * 2143 ifelse(q >= location, kappa, 1/kappa) 2144 temp5 <- exp(exponent) / (1 + kappa^2) 2145 index1 <- (q < location) 2146 2147 2148 if (lower.tail) { 2149 if (log.p) { 2150 ans <- log1p(-exp(exponent) / (1 + kappa^2)) 2151 logtemp5 <- exponent - log1p(kappa^2) 2152 ans[index1] <- 2 * log(kappa[index1]) + logtemp5[index1] 2153 } else { 2154 ans <- (kappa^2 - expm1(exponent)) / (1 + kappa^2) 2155 ans[index1] <- (kappa[index1])^2 * temp5[index1] 2156 } 2157 } else { 2158 if (log.p) { 2159 ans <- exponent - log1p(kappa^2) # logtemp5 2160 ans[index1] <- log1p(-(kappa[index1])^2 * temp5[index1]) 2161 } else { 2162 ans <- temp5 2163 ans[index1] <- (1 + (kappa[index1])^2 * 2164 (-expm1(exponent[index1]))) / (1+(kappa[index1])^2) 2165 } 2166 } 2167 indexTF <- (scale > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2168 ans[!indexTF] <- NaN 2169 ans 2170} 2171 2172 2173 2174qalap <- function(p, location = 0, scale = 1, tau = 0.5, 2175 kappa = sqrt(tau / (1 - tau)), 2176 lower.tail = TRUE, log.p = FALSE) { 2177 2178 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 2179 stop("bad input for argument 'lower.tail'") 2180 2181 if (!is.logical(log.p) || length(log.p) != 1) 2182 stop("bad input for argument 'log.p'") 2183 2184 NN <- max(length(p), length(location), length(scale), length(kappa), 2185 length(tau)) 2186 if (length(p) != NN) p <- rep_len(p, NN) 2187 if (length(location) != NN) location <- rep_len(location, NN) 2188 if (length(scale) != NN) scale <- rep_len(scale, NN) 2189 if (length(kappa) != NN) kappa <- rep_len(kappa, NN) 2190 if (length(tau) != NN) tau <- rep_len(tau, NN) 2191 2192 2193 2194 temp5 <- kappa^2 / (1 + kappa^2) 2195 if (lower.tail) { 2196 if (log.p) { 2197 ans <- exp(p) 2198 index1 <- (exp(p) <= temp5) 2199 exponent <- exp(p[index1]) / temp5[index1] 2200 ans[index1] <- location[index1] + (scale[index1] * kappa[index1]) * 2201 log(exponent) / sqrt(2) 2202 ans[!index1] <- location[!index1] - 2203 (scale[!index1] / kappa[!index1]) * 2204 (log1p((kappa[!index1])^2) + 2205 log(-expm1(p[!index1]))) / sqrt(2) 2206 } else { 2207 ans <- p 2208 index1 <- (p <= temp5) 2209 exponent <- p[index1] / temp5[index1] 2210 ans[index1] <- location[index1] + (scale[index1] * kappa[index1]) * 2211 log(exponent) / sqrt(2) 2212 ans[!index1] <- location[!index1] - 2213 (scale[!index1] / kappa[!index1]) * 2214 (log1p((kappa[!index1])^2) + 2215 log1p(-p[!index1])) / sqrt(2) 2216 } 2217 } else { 2218 if (log.p) { 2219 ans <- -expm1(p) 2220 index1 <- (-expm1(p) <= temp5) 2221 exponent <- -expm1(p[index1]) / temp5[index1] 2222 ans[index1] <- location[index1] + (scale[index1] * kappa[index1]) * 2223 log(exponent) / sqrt(2) 2224 ans[!index1] <- location[!index1] - 2225 (scale[!index1] / kappa[!index1]) * 2226 (log1p((kappa[!index1])^2) + 2227 p[!index1]) / sqrt(2) 2228 } else { 2229 ans <- exp(log1p(-p)) 2230 index1 <- (p >= (1 / (1+kappa^2))) 2231 exponent <- exp(log1p(-p[index1])) / temp5[index1] 2232 ans[index1] <- location[index1] + 2233 (scale[index1] * kappa[index1]) * log(exponent) / sqrt(2) 2234 ans[!index1] <- location[!index1] - 2235 (scale[!index1] / kappa[!index1]) * 2236 (log1p((kappa[!index1])^2) + 2237 log(p[!index1])) / sqrt(2) 2238 } 2239 } 2240 2241 indexTF <- (scale > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2242 ans[!indexTF] <- NaN 2243 ans 2244} 2245 2246 2247 2248 2249 2250 2251 2252rloglap <- function(n, location.ald = 0, scale.ald = 1, tau = 0.5, 2253 kappa = sqrt(tau/(1-tau))) { 2254 use.n <- if ((length.n <- length(n)) > 1) length.n else 2255 if (!is.Numeric(n, integer.valued = TRUE, 2256 length.arg = 1, positive = TRUE)) 2257 stop("bad input for argument 'n'") else n 2258 location.ald <- rep_len(location.ald, use.n) 2259 scale.ald <- rep_len(scale.ald, use.n) 2260 tau <- rep_len(tau, use.n) 2261 kappa <- rep_len(kappa, use.n) 2262 ans <- exp(location.ald) * 2263 (runif(use.n)^kappa / runif(use.n)^(1/kappa))^(scale.ald / sqrt(2)) 2264 indexTF <- (scale.ald > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2265 ans[!indexTF] <- NaN 2266 ans 2267} 2268 2269 2270 2271dloglap <- function(x, location.ald = 0, scale.ald = 1, tau = 0.5, 2272 kappa = sqrt(tau/(1-tau)), log = FALSE) { 2273 if (!is.logical(log.arg <- log) || length(log) != 1) 2274 stop("bad input for argument 'log'") 2275 rm(log) 2276 2277 2278 scale <- scale.ald 2279 location <- location.ald 2280 NN <- max(length(x), length(location), 2281 length(scale), length(kappa), length(tau)) 2282 2283 if (length(x) != NN) x <- rep_len(x, NN) 2284 if (length(location) != NN) location <- rep_len(location, NN) 2285 if (length(scale) != NN) scale <- rep_len(scale, NN) 2286 if (length(kappa) != NN) kappa <- rep_len(kappa, NN) 2287 if (length(tau) != NN) tau <- rep_len(tau, NN) 2288 2289 2290 Alpha <- sqrt(2) * kappa / scale.ald 2291 Beta <- sqrt(2) / (scale.ald * kappa) 2292 Delta <- exp(location.ald) 2293 exponent <- ifelse(x >= Delta, -(Alpha+1), (Beta-1)) * 2294 (log(x) - location.ald) 2295 logdensity <- -location.ald + log(Alpha) + log(Beta) - 2296 log(Alpha + Beta) + exponent 2297 indexTF <- (scale.ald > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2298 logdensity[!indexTF] <- NaN 2299 logdensity[x < 0 & indexTF] <- -Inf 2300 if (log.arg) logdensity else exp(logdensity) 2301} 2302 2303 2304 2305qloglap <- function(p, location.ald = 0, scale.ald = 1, 2306 tau = 0.5, kappa = sqrt(tau/(1-tau)), 2307 lower.tail = TRUE, log.p = FALSE) { 2308 2309 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 2310 stop("bad input for argument 'lower.tail'") 2311 2312 if (!is.logical(log.p) || length(log.p) != 1) 2313 stop("bad input for argument 'log.p'") 2314 2315 2316 NN <- max(length(p), length(location.ald), length(scale.ald), 2317 length(kappa)) 2318 p <- rep_len(p, NN) 2319 location <- rep_len(location.ald, NN) 2320 scale <- rep_len(scale.ald, NN) 2321 kappa <- rep_len(kappa, NN) 2322 tau <- rep_len(tau, NN) 2323 2324 2325 Alpha <- sqrt(2) * kappa / scale.ald 2326 Beta <- sqrt(2) / (scale.ald * kappa) 2327 Delta <- exp(location.ald) 2328 temp9 <- Alpha + Beta 2329 2330 2331 if (lower.tail) { 2332 if (log.p) { 2333 ln.p <- p 2334 ans <- ifelse((exp(ln.p) > Alpha / temp9), 2335 Delta * (-expm1(ln.p) * temp9 / Beta)^(-1/Alpha), 2336 Delta * (exp(ln.p) * temp9 / Alpha)^(1/Beta)) 2337 ans[ln.p > 0] <- NaN 2338 } else { 2339 ans <- ifelse((p > Alpha / temp9), 2340 Delta * exp((-1/Alpha) * (log1p(-p) + 2341 log(temp9/Beta))), 2342 Delta * (p * temp9 / Alpha)^(1/Beta)) 2343 ans[p < 0] <- NaN 2344 ans[p == 0] <- 0 2345 ans[p == 1] <- Inf 2346 ans[p > 1] <- NaN 2347 } 2348 } else { 2349 if (log.p) { 2350 ln.p <- p 2351 ans <- ifelse((-expm1(ln.p) > Alpha / temp9), 2352 Delta * (exp(ln.p) * temp9 / Beta)^(-1/Alpha), 2353 Delta * (-expm1(ln.p) * temp9 / Alpha)^(1/Beta)) 2354 ans[ln.p > 0] <- NaN 2355 } else { 2356 ans <- ifelse((p < (temp9 - Alpha) / temp9), 2357 Delta * (p * temp9 / Beta)^(-1/Alpha), 2358 Delta * exp((1/Beta)*(log1p(-p) + log(temp9/Alpha)))) 2359 ans[p < 0] <- NaN 2360 ans[p == 0] <- Inf 2361 ans[p == 1] <- 0 2362 ans[p > 1] <- NaN 2363 } 2364 } 2365 indexTF <- (scale.ald > 0) & (tau > 0) & (tau < 1) & (kappa > 0) 2366 ans[!indexTF] <- NaN 2367 ans 2368} 2369 2370 2371 2372ploglap <- function(q, location.ald = 0, scale.ald = 1, 2373 tau = 0.5, kappa = sqrt(tau/(1-tau)), 2374 lower.tail = TRUE, log.p = FALSE) { 2375 2376 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 2377 stop("bad input for argument 'lower.tail'") 2378 2379 if (!is.logical(log.p) || length(log.p) != 1) 2380 stop("bad input for argument 'log.p'") 2381 2382 NN <- max(length(q), length(location.ald), length(scale.ald), 2383 length(kappa)) 2384 location <- rep_len(location.ald, NN) 2385 scale <- rep_len(scale.ald, NN) 2386 kappa <- rep_len(kappa, NN) 2387 q <- rep_len(q, NN) 2388 tau <- rep_len(tau, NN) 2389 2390 Alpha <- sqrt(2) * kappa / scale.ald 2391 Beta <- sqrt(2) / (scale.ald * kappa) 2392 Delta <- exp(location.ald) 2393 2394 temp9 <- Alpha + Beta 2395 index1 <- (Delta <= q) 2396 2397 2398 if (lower.tail) { 2399 if (log.p) { 2400 ans <- log((Alpha / temp9) * (q / Delta)^(Beta)) 2401 ans[index1] <- log1p((-(Beta/temp9) * (Delta/q)^(Alpha))[index1]) 2402 ans[q <= 0 ] <- -Inf 2403 ans[q == Inf] <- 0 2404 } else { 2405 ans <- (Alpha / temp9) * (q / Delta)^(Beta) 2406 ans[index1] <- -expm1((log(Beta/temp9) + 2407 Alpha * log(Delta/q)))[index1] 2408 ans[q <= 0] <- 0 2409 ans[q == Inf] <- 1 2410 } 2411 } else { 2412 if (log.p) { 2413 ans <- log1p(-(Alpha / temp9) * (q / Delta)^(Beta)) 2414 ans[index1] <- log(((Beta/temp9) * (Delta/q)^(Alpha))[index1]) 2415 ans[q <= 0] <- 0 2416 ans[q == Inf] <- -Inf 2417 } else { 2418 ans <- -expm1(log(Alpha/temp9) + Beta * log(q/Delta)) 2419 ans[index1] <- ((Beta/temp9) * (Delta/q)^(Alpha))[index1] 2420 ans[q <= 0] <- 1 2421 ans[q == Inf] <- 0 2422 } 2423 } 2424 2425 indexTF <- (scale.ald > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2426 ans[!indexTF] <- NaN 2427 ans 2428} 2429 2430 2431 2432 2433rlogitlap <- function(n, location.ald = 0, scale.ald = 1, tau = 0.5, 2434 kappa = sqrt(tau/(1-tau))) { 2435 logitlink(ralap(n = n, location = location.ald, scale = scale.ald, 2436 tau = tau, kappa = kappa), 2437 inverse = TRUE) # earg = earg 2438} 2439 2440 2441 2442dlogitlap <- function(x, location.ald = 0, scale.ald = 1, tau = 0.5, 2443 kappa = sqrt(tau/(1-tau)), log = FALSE) { 2444 if (!is.logical(log.arg <- log) || length(log) != 1) 2445 stop("bad input for argument 'log'") 2446 rm(log) 2447 2448 2449 2450 NN <- max(length(x), length(location.ald), 2451 length(scale.ald), length(kappa)) 2452 location <- rep_len(location.ald, NN) 2453 scale <- rep_len(scale.ald, NN) 2454 kappa <- rep_len(kappa, NN) 2455 x <- rep_len(x, NN) 2456 tau <- rep_len(tau, NN) 2457 2458 Alpha <- sqrt(2) * kappa / scale.ald 2459 Beta <- sqrt(2) / (scale.ald * kappa) 2460 Delta <- logitlink(location.ald, inverse = TRUE) # earg = earg 2461 2462 exponent <- ifelse(x >= Delta, -Alpha, Beta) * 2463 (logitlink(x) - # earg = earg 2464 location.ald) 2465 logdensity <- log(Alpha) + log(Beta) - log(Alpha + Beta) - 2466 log(x) - log1p(-x) + exponent 2467 indexTF <- (scale.ald > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2468 logdensity[!indexTF] <- NaN 2469 logdensity[x < 0 & indexTF] <- -Inf 2470 logdensity[x > 1 & indexTF] <- -Inf 2471 if (log.arg) logdensity else exp(logdensity) 2472} 2473 2474 2475 2476qlogitlap <- function(p, location.ald = 0, scale.ald = 1, 2477 tau = 0.5, kappa = sqrt(tau/(1-tau))) { 2478 qqq <- qalap(p = p, location = location.ald, scale = scale.ald, 2479 tau = tau, kappa = kappa) 2480 ans <- logitlink(qqq, inverse = TRUE) # earg = earg 2481 ans[(p < 0) | (p > 1)] <- NaN 2482 ans[p == 0] <- 0 2483 ans[p == 1] <- 1 2484 ans 2485} 2486 2487 2488 2489plogitlap <- function(q, location.ald = 0, scale.ald = 1, 2490 tau = 0.5, kappa = sqrt(tau/(1-tau))) { 2491 NN <- max(length(q), length(location.ald), length(scale.ald), 2492 length(kappa)) 2493 location.ald <- rep_len(location.ald, NN) 2494 scale.ald <- rep_len(scale.ald, NN) 2495 kappa <- rep_len(kappa, NN) 2496 q <- rep_len(q, NN) 2497 tau <- rep_len(tau, NN) 2498 2499 indexTF <- (q > 0) & (q < 1) 2500 qqq <- logitlink(q[indexTF]) # earg = earg 2501 ans <- q 2502 ans[indexTF] <- palap(q = qqq, location = location.ald[indexTF], 2503 scale = scale.ald[indexTF], 2504 tau = tau[indexTF], kappa = kappa[indexTF]) 2505 ans[q >= 1] <- 1 2506 ans[q <= 0] <- 0 2507 ans 2508} 2509 2510 2511 2512 2513rprobitlap <- function(n, location.ald = 0, scale.ald = 1, tau = 0.5, 2514 kappa = sqrt(tau/(1-tau))) { 2515 2516 2517 2518 probitlink(ralap(n = n, location = location.ald, scale = scale.ald, 2519 tau = tau, kappa = kappa), 2520 inverse = TRUE) 2521} 2522 2523 2524 2525dprobitlap <- 2526 function(x, location.ald = 0, scale.ald = 1, tau = 0.5, 2527 kappa = sqrt(tau/(1-tau)), log = FALSE, 2528 meth2 = TRUE) { 2529 if (!is.logical(log.arg <- log) || length(log) != 1) 2530 stop("bad input for argument 'log'") 2531 rm(log) 2532 2533 2534 2535 NN <- max(length(x), length(location.ald), length(scale.ald), 2536 length(kappa)) 2537 location.ald <- rep_len(location.ald, NN) 2538 scale.ald <- rep_len(scale.ald, NN) 2539 kappa <- rep_len(kappa, NN) 2540 x <- rep_len(x, NN) 2541 tau <- rep_len(tau, NN) 2542 2543 logdensity <- x * NaN 2544 index1 <- (x > 0) & (x < 1) 2545 indexTF <- (scale.ald > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2546 if (meth2) { 2547 dx.dy <- x 2548 use.x <- probitlink(x[index1]) # earg = earg 2549 logdensity[index1] <- 2550 dalap(x = use.x, location = location.ald[index1], 2551 scale = scale.ald[index1], tau = tau[index1], 2552 kappa = kappa[index1], log = TRUE) 2553 } else { 2554 Alpha <- sqrt(2) * kappa / scale.ald 2555 Beta <- sqrt(2) / (scale.ald * kappa) 2556 Delta <- pnorm(location.ald) 2557 use.x <- qnorm(x) # qnorm(x[index1]) 2558 log.dy.dw <- dnorm(use.x, log = TRUE) 2559 2560 exponent <- ifelse(x >= Delta, -Alpha, Beta) * 2561 (use.x - location.ald) - log.dy.dw 2562 2563 logdensity[index1] <- (log(Alpha) + log(Beta) - 2564 log(Alpha + Beta) + exponent)[index1] 2565 } 2566 logdensity[!indexTF] <- NaN 2567 logdensity[x < 0 & indexTF] <- -Inf 2568 logdensity[x > 1 & indexTF] <- -Inf 2569 2570 if (meth2) { 2571 dx.dy[index1] <- probitlink(x[index1], # earg = earg, 2572 inverse = TRUE, 2573 deriv = 1) 2574 dx.dy[!index1] <- 0 2575 dx.dy[!indexTF] <- NaN 2576 if (log.arg) logdensity - log(abs(dx.dy)) else 2577 exp(logdensity) / abs(dx.dy) 2578 } else { 2579 if (log.arg) logdensity else exp(logdensity) 2580 } 2581} 2582 2583 2584qprobitlap <- function(p, location.ald = 0, scale.ald = 1, 2585 tau = 0.5, kappa = sqrt(tau/(1-tau))) { 2586 qqq <- qalap(p = p, location = location.ald, scale = scale.ald, 2587 tau = tau, kappa = kappa) 2588 ans <- probitlink(qqq, inverse = TRUE) # , earg = earg 2589 ans[(p < 0) | (p > 1)] = NaN 2590 ans[p == 0] <- 0 2591 ans[p == 1] <- 1 2592 ans 2593} 2594 2595 2596 2597pprobitlap <- function(q, location.ald = 0, scale.ald = 1, 2598 tau = 0.5, kappa = sqrt(tau/(1-tau))) { 2599 NN <- max(length(q), length(location.ald), length(scale.ald), 2600 length(kappa)) 2601 location.ald <- rep_len(location.ald, NN) 2602 scale.ald <- rep_len(scale.ald, NN) 2603 kappa <- rep_len(kappa, NN) 2604 q <- rep_len(q, NN) 2605 tau <- rep_len(tau, NN) 2606 2607 indexTF <- (q > 0) & (q < 1) 2608 qqq <- probitlink(q[indexTF]) # earg = earg 2609 ans <- q 2610 ans[indexTF] <- palap(q = qqq, location = location.ald[indexTF], 2611 scale = scale.ald[indexTF], 2612 tau = tau[indexTF], kappa = kappa[indexTF]) 2613 ans[q >= 1] <- 1 2614 ans[q <= 0] <- 0 2615 ans 2616} 2617 2618 2619 2620 2621 2622rclogloglap <- function(n, location.ald = 0, scale.ald = 1, tau = 0.5, 2623 kappa = sqrt(tau/(1-tau))) { 2624 clogloglink(ralap(n = n, location = location.ald, scale = scale.ald, 2625 tau = tau, kappa = kappa), # earg = earg, 2626 inverse = TRUE) 2627} 2628 2629 2630 2631dclogloglap <- function(x, location.ald = 0, scale.ald = 1, tau = 0.5, 2632 kappa = sqrt(tau/(1-tau)), log = FALSE, 2633 meth2 = TRUE) { 2634 if (!is.logical(log.arg <- log) || length(log) != 1) 2635 stop("bad input for argument 'log'") 2636 rm(log) 2637 2638 2639 2640 NN <- max(length(x), length(location.ald), length(scale.ald), 2641 length(kappa)) 2642 location.ald <- rep_len(location.ald, NN) 2643 scale.ald <- rep_len(scale.ald, NN) 2644 kappa <- rep_len(kappa, NN) 2645 x <- rep_len(x, NN) 2646 tau <- rep_len(tau, NN) 2647 2648 logdensity <- x * NaN 2649 index1 <- (x > 0) & (x < 1) 2650 indexTF <- (scale.ald > 0) & (tau > 0) & (tau < 1) & (kappa > 0) # & 2651 if (meth2) { 2652 dx.dy <- x 2653 use.w <- clogloglink(x[index1]) # earg = earg 2654 logdensity[index1] <- 2655 dalap(x = use.w, location = location.ald[index1], 2656 scale = scale.ald[index1], 2657 tau = tau[index1], 2658 kappa = kappa[index1], log = TRUE) 2659 2660 } else { 2661 Alpha <- sqrt(2) * kappa / scale.ald 2662 Beta <- sqrt(2) / (scale.ald * kappa) 2663 Delta <- clogloglink(location.ald, inverse = TRUE) 2664 2665 exponent <- ifelse(x >= Delta, -(Alpha+1), Beta-1) * 2666 log(-log1p(-x)) + 2667 ifelse(x >= Delta, Alpha, -Beta) * location.ald 2668 logdensity[index1] <- (log(Alpha) + log(Beta) - 2669 log(Alpha + Beta) - log1p(-x) + exponent)[index1] 2670 } 2671 logdensity[!indexTF] <- NaN 2672 logdensity[x < 0 & indexTF] <- -Inf 2673 logdensity[x > 1 & indexTF] <- -Inf 2674 2675 if (meth2) { 2676 dx.dy[index1] <- clogloglink(x[index1], # earg = earg, 2677 inverse = TRUE, deriv = 1) 2678 dx.dy[!index1] <- 0 2679 dx.dy[!indexTF] <- NaN 2680 if (log.arg) logdensity - log(abs(dx.dy)) else 2681 exp(logdensity) / abs(dx.dy) 2682 } else { 2683 if (log.arg) logdensity else exp(logdensity) 2684 } 2685} 2686 2687 2688 2689qclogloglap <- function(p, location.ald = 0, scale.ald = 1, 2690 tau = 0.5, kappa = sqrt(tau/(1-tau))) { 2691 qqq <- qalap(p = p, location = location.ald, scale = scale.ald, 2692 tau = tau, kappa = kappa) 2693 ans <- clogloglink(qqq, inverse = TRUE) # , earg = earg 2694 ans[(p < 0) | (p > 1)] <- NaN 2695 ans[p == 0] <- 0 2696 ans[p == 1] <- 1 2697 ans 2698} 2699 2700 2701 2702pclogloglap <- function(q, location.ald = 0, scale.ald = 1, 2703 tau = 0.5, kappa = sqrt(tau/(1-tau))) { 2704 NN <- max(length(q), length(location.ald), length(scale.ald), 2705 length(kappa)) 2706 location.ald <- rep_len(location.ald, NN) 2707 scale.ald <- rep_len(scale.ald, NN) 2708 kappa <- rep_len(kappa, NN) 2709 q <- rep_len(q, NN) 2710 tau <- rep_len(tau, NN) 2711 2712 2713 indexTF <- (q > 0) & (q < 1) 2714 qqq <- clogloglink(q[indexTF]) # earg = earg 2715 ans <- q 2716 ans[indexTF] <- palap(q = qqq, location = location.ald[indexTF], 2717 scale = scale.ald[indexTF], 2718 tau = tau[indexTF], kappa = kappa[indexTF]) 2719 ans[q >= 1] <- 1 2720 ans[q <= 0] <- 0 2721 ans 2722} 2723 2724 2725 2726 2727 2728 2729 2730 2731 2732 2733 2734 2735alaplace2.control <- function(maxit = 100, ...) { 2736 list(maxit = maxit) 2737} 2738 2739 2740 alaplace2 <- 2741 function(tau = NULL, 2742 llocation = "identitylink", lscale = "loglink", 2743 ilocation = NULL, iscale = NULL, 2744 kappa = sqrt(tau / (1-tau)), 2745 ishrinkage = 0.95, 2746 2747 parallel.locat = TRUE ~ 0, 2748 2749 parallel.scale = FALSE ~ 0, 2750 2751 digt = 4, 2752 idf.mu = 3, 2753 imethod = 1, 2754 zero = "scale") { 2755 2756 2757 2758 apply.parint.locat <- FALSE 2759 apply.parint.scale <- TRUE 2760 2761 2762 2763 2764 llocat <- as.list(substitute(llocation)) 2765 elocat <- link2list(llocat) 2766 llocat <- attr(elocat, "function.name") 2767 2768 lscale <- as.list(substitute(lscale)) 2769 escale <- link2list(lscale) 2770 lscale <- attr(escale, "function.name") 2771 2772 ilocat <- ilocation 2773 2774 2775 2776 if (!is.Numeric(kappa, positive = TRUE)) 2777 stop("bad input for argument 'kappa'") 2778 if (!is.Numeric(imethod, length.arg = 1, 2779 integer.valued = TRUE, positive = TRUE) || 2780 imethod > 4) 2781 stop("argument 'imethod' must be 1, 2 or ... 4") 2782 if (length(iscale) && 2783 !is.Numeric(iscale, positive = TRUE)) 2784 stop("bad input for argument 'iscale'") 2785 if (!is.Numeric(ishrinkage, length.arg = 1) || 2786 ishrinkage < 0 || 2787 ishrinkage > 1) 2788 stop("bad input for argument 'ishrinkage'") 2789 if (length(tau) && 2790 max(abs(kappa - sqrt(tau / (1 - tau)))) > 1.0e-6) 2791 stop("arguments 'kappa' and 'tau' do not match") 2792 2793 2794 2795 fittedMean <- FALSE 2796 if (!is.logical(fittedMean) || length(fittedMean) != 1) 2797 stop("bad input for argument 'fittedMean'") 2798 2799 2800 2801 new("vglmff", 2802 blurb = c("Two-parameter asymmetric Laplace distribution\n\n", 2803 "Links: ", 2804 namesof("location", llocat, earg = elocat), ", ", 2805 namesof("scale", lscale, earg = escale), # ", ", 2806 "\n\n", 2807 "Mean: ", 2808 "location + scale * (1/kappa - kappa) / sqrt(2)", "\n", 2809 "Quantiles: location", "\n", 2810 "Variance: scale^2 * (1 + kappa^4) / (2 * kappa^2)"), 2811 2812 2813 2814 2815 constraints = eval(substitute(expression({ 2816 2817 2818 onemat <- matrix(1, Mdiv2, 1) 2819 constraints.orig <- constraints 2820 2821 2822 cm1.locat <- kronecker(diag(Mdiv2), rbind(1, 0)) 2823 cmk.locat <- kronecker(onemat, rbind(1, 0)) 2824 con.locat <- cm.VGAM(cmk.locat, 2825 x = x, bool = .parallel.locat , 2826 constraints = constraints.orig, 2827 apply.int = .apply.parint.locat , 2828 cm.default = cm1.locat, 2829 cm.intercept.default = cm1.locat) 2830 2831 2832 2833 cm1.scale <- kronecker(diag(Mdiv2), rbind(0, 1)) 2834 cmk.scale <- kronecker(onemat, rbind(0, 1)) 2835 con.scale <- cm.VGAM(cmk.scale, 2836 x = x, bool = .parallel.scale , 2837 constraints = constraints.orig, 2838 apply.int = .apply.parint.scale , 2839 cm.default = cm1.scale, 2840 cm.intercept.default = cm1.scale) 2841 2842 con.use <- con.scale 2843 for (klocal in seq_along(con.scale)) { 2844 con.use[[klocal]] <- cbind(con.locat[[klocal]], 2845 con.scale[[klocal]]) 2846 } 2847 2848 2849 constraints <- con.use 2850 2851 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 2852 predictors.names = predictors.names, 2853 M1 = M1) 2854 }), list( .parallel.locat = parallel.locat, 2855 .parallel.scale = parallel.scale, 2856 .zero = zero, 2857 .apply.parint.scale = apply.parint.scale, 2858 .apply.parint.locat = apply.parint.locat ))), 2859 2860 2861 infos = eval(substitute(function(...) { 2862 list(M1 = 2, 2863 Q1 = 1, 2864 summary.pvalues = FALSE, 2865 expected = TRUE, # 20161117 2866 multipleResponses = TRUE, # FALSE, 2867 parameters.names = c("location", "scale"), 2868 true.mu = .fittedMean , 2869 zero = .zero , 2870 tau = .tau , 2871 kappa = .kappa ) 2872 }, list( .tau = tau, 2873 .kappa = kappa, 2874 .fittedMean = fittedMean, 2875 .zero = zero ))), 2876 2877 2878 initialize = eval(substitute(expression({ 2879 M1 <- 2 2880 2881 temp5 <- 2882 w.y.check(w = w, y = y, 2883 ncol.w.max = if (length( .kappa ) > 1) 1 else Inf, 2884 ncol.y.max = if (length( .kappa ) > 1) 1 else Inf, 2885 out.wy = TRUE, 2886 colsyperw = 1, # Uncommented out 20140621 2887 maximize = TRUE) 2888 w <- temp5$w 2889 y <- temp5$y 2890 2891 2892 extra$ncoly <- ncoly <- ncol(y) 2893 if ((ncoly > 1) && (length( .kappa ) > 1)) 2894 stop("response must be a vector if 'kappa' or 'tau' ", 2895 "has a length greater than one") 2896 2897 2898 2899 extra$kappa <- .kappa 2900 extra$tau <- extra$kappa^2 / (1 + extra$kappa^2) 2901 2902 extra$Mdiv2 <- Mdiv2 <- max(ncoly, length( .kappa )) 2903 extra$M <- M <- M1 * Mdiv2 2904 extra$n <- n 2905 2906 2907 2908 extra$tau.names <- tau.names <- 2909 paste("(tau = ", round(extra$tau, digits = .digt), ")", sep = "") 2910 extra$Y.names <- Y.names <- if (ncoly > 1) dimnames(y)[[2]] else "y" 2911 if (is.null(Y.names) || any(Y.names == "")) 2912 extra$Y.names <- Y.names <- paste("y", 1:ncoly, sep = "") 2913 extra$y.names <- y.names <- 2914 if (ncoly > 1) paste(Y.names, tau.names, sep = "") else tau.names 2915 2916 extra$individual <- FALSE 2917 2918 2919 mynames1 <- param.names("location", Mdiv2, skip1 = TRUE) 2920 mynames2 <- param.names("scale", Mdiv2, skip1 = TRUE) 2921 predictors.names <- 2922 c(namesof(mynames1, .llocat , earg = .elocat , tag = FALSE), 2923 namesof(mynames2, .lscale , earg = .escale , tag = FALSE)) 2924 predictors.names <- 2925 predictors.names[interleave.VGAM(M, M1 = M1)] 2926 2927 2928 2929 2930 locat.init <- scale.init <- matrix(0, n, Mdiv2) 2931 if (!length(etastart)) { 2932 for (jay in 1:Mdiv2) { 2933 y.use <- if (ncoly > 1) y[, jay] else y 2934 Jay <- if (ncoly > 1) jay else 1 2935 if ( .imethod == 1) { 2936 locat.init[, jay] <- weighted.mean(y.use, w[, Jay]) 2937 scale.init[, jay] <- sqrt(var(y.use) / 2) 2938 } else if ( .imethod == 2) { 2939 locat.init[, jay] <- median(y.use) 2940 scale.init[, jay] <- sqrt(sum(c(w[, Jay]) * 2941 abs(y - median(y.use))) / (sum(w[, Jay]) * 2)) 2942 } else if ( .imethod == 3) { 2943 Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 2944 y = y.use, w = w[, Jay], 2945 df = .idf.mu ) 2946 locat.init[, jay] <- predict(Fit5, x = x[, min(ncol(x), 2)])$y 2947 scale.init[, jay] <- sqrt(sum(c(w[, Jay]) * 2948 abs(y.use - median(y.use))) / ( 2949 sum(w[, Jay]) * 2)) 2950 } else { 2951 use.this <- weighted.mean(y.use, w[, Jay]) 2952 locat.init[, jay] <- (1 - .ishrinkage ) * y.use + 2953 .ishrinkage * use.this 2954 scale.init[, jay] <- 2955 sqrt(sum(c(w[, Jay]) * 2956 abs(y.use - median(y.use ))) / (sum(w[, Jay]) * 2)) 2957 } 2958 } 2959 2960 2961 2962 if (length( .ilocat )) { 2963 locat.init <- matrix( .ilocat , n, Mdiv2, byrow = TRUE) 2964 } 2965 if (length( .iscale )) { 2966 scale.init <- matrix( .iscale , n, Mdiv2, byrow = TRUE) 2967 } 2968 2969 etastart <- 2970 cbind(theta2eta(locat.init, .llocat , earg = .elocat ), 2971 theta2eta(scale.init, .lscale , earg = .escale )) 2972 etastart <- etastart[, interleave.VGAM(M, M1 = M1), drop = FALSE] 2973 } 2974 }), list( .imethod = imethod, 2975 .idf.mu = idf.mu, 2976 .ishrinkage = ishrinkage, .digt = digt, 2977 .elocat = elocat, .escale = escale, 2978 .llocat = llocat, .lscale = lscale, .kappa = kappa, 2979 .ilocat = ilocat, .iscale = iscale ))), 2980 2981 linkinv = eval(substitute(function(eta, extra = NULL) { 2982 M1 <- 2 2983 Mdiv2 <- ncol(eta) / M1 # extra$Mdiv2 2984 vTF <- c(TRUE, FALSE) 2985 locat <- eta2theta(eta[, vTF, drop = FALSE], .llocat , 2986 earg = .elocat ) 2987 dimnames(locat) <- list(dimnames(eta)[[1]], extra$y.names) 2988 myans <- if ( .fittedMean ) { 2989 kappamat <- matrix(extra$kappa, extra$n, extra$Mdiv2, 2990 byrow = TRUE) 2991 Scale <- eta2theta(eta[, !vTF, drop = FALSE], .lscale , 2992 earg = .escale ) 2993 locat + Scale * (1/kappamat - kappamat) 2994 } else { 2995 locat 2996 } 2997 dimnames(myans) <- list(dimnames(myans)[[1]], extra$y.names) 2998 myans 2999 }, list( .elocat = elocat, .llocat = llocat, 3000 .escale = escale, .lscale = lscale, 3001 .fittedMean = fittedMean, 3002 .kappa = kappa ))), 3003 last = eval(substitute(expression({ 3004 M1 <- 2 # extra$M1 3005 Mdiv2 <- ncol(eta) / M1 # extra$Mdiv2 3006 3007 3008 3009 misc$link <- setNames(c(rep_len( .llocat , Mdiv2), 3010 rep_len( .lscale , Mdiv2)), 3011 c(mynames1, mynames2))[interleave.VGAM(M, M1 = M1)] 3012 3013 3014 3015 misc$earg <- vector("list", M) 3016 for (ii in 1:Mdiv2) { 3017 misc$earg[[M1 * ii - 1]] <- .elocat 3018 misc$earg[[M1 * ii ]] <- .escale 3019 } 3020 names(misc$earg) <- names(misc$link) 3021 3022 3023 extra$kappa <- misc$kappa <- .kappa 3024 extra$tau <- misc$tau <- misc$kappa^2 / (1 + misc$kappa^2) 3025 3026 extra$percentile <- numeric(Mdiv2) # length(misc$kappa) 3027 locat <- as.matrix(locat) 3028 for (ii in 1:Mdiv2) { 3029 y.use <- if (ncoly > 1) y[, ii] else y 3030 Jay <- if (ncoly > 1) ii else 1 3031 extra$percentile[ii] <- 100 * weighted.mean(y.use <= locat[, ii], 3032 w[, Jay]) 3033 } 3034 names(extra$percentile) <- y.names 3035 }), list( .elocat = elocat, .llocat = llocat, 3036 .escale = escale, .lscale = lscale, 3037 .fittedMean = fittedMean, 3038 .kappa = kappa ))), 3039 3040 loglikelihood = eval(substitute( 3041 function(mu, y, w, residuals = FALSE, eta, 3042 extra = NULL, 3043 summation = TRUE) { 3044 M1 <- 2 3045 Mdiv2 <- ncol(eta) / M1 # extra$Mdiv2 3046 ymat <- matrix(y, extra$n, extra$Mdiv2) 3047 kappamat <- matrix(extra$kappa, extra$n, extra$Mdiv2, byrow = TRUE) 3048 3049 vTF <- c(TRUE, FALSE) 3050 locat <- eta2theta(eta[, vTF, drop = FALSE], .llocat , 3051 earg = .elocat ) 3052 Scale <- eta2theta(eta[, !vTF, drop = FALSE], .lscale , 3053 earg = .escale ) 3054 if (residuals) { 3055 stop("loglikelihood residuals not implemented yet") 3056 } else { 3057 ll.elts <- c(w) * dalap(x = c(ymat), location = c(locat), 3058 scale = c(Scale), kappa = c(kappamat), 3059 log = TRUE) 3060 if (summation) { 3061 sum(ll.elts) 3062 } else { 3063 ll.elts 3064 } 3065 } 3066 }, list( .elocat = elocat, .llocat = llocat, 3067 .escale = escale, .lscale = lscale, 3068 .kappa = kappa ))), 3069 vfamily = c("alaplace2"), 3070 validparams = eval(substitute(function(eta, y, extra = NULL) { 3071 vTF <- c(TRUE, FALSE) 3072 locat <- eta2theta(eta[, vTF, drop = FALSE], .llocat , 3073 earg = .elocat ) 3074 Scale <- eta2theta(eta[, !vTF, drop = FALSE], .lscale , 3075 earg = .escale ) 3076 okay1 <- all(is.finite(locat)) && 3077 all(is.finite(Scale)) && all(0 < Scale) 3078 okay1 3079 }, list( .elocat = elocat, .llocat = llocat, 3080 .escale = escale, .lscale = lscale, 3081 .kappa = kappa ))), 3082 3083 3084 3085 3086 simslot = eval(substitute( 3087 function(object, nsim) { 3088 3089 pwts <- if (length(pwts <- object@prior.weights) > 0) 3090 pwts else weights(object, type = "prior") 3091 if (any(pwts != 1)) 3092 warning("ignoring prior weights") 3093 eta <- predict(object) 3094 extra <- object@extra 3095 vTF <- c(TRUE, FALSE) 3096 locat <- eta2theta(eta[, vTF, drop = FALSE], .llocat , earg = .elocat ) 3097 Scale <- eta2theta(eta[, !vTF, drop = FALSE], .lscale , earg = .escale ) 3098 kappamat <- matrix(extra$kappa, extra$n, extra$Mdiv2, byrow = TRUE) 3099 ralap(nsim * length(Scale), location = c(locat), 3100 scale = c(Scale), kappa = c(kappamat)) 3101 }, list( .elocat = elocat, .llocat = llocat, 3102 .escale = escale, .lscale = lscale, 3103 .kappa = kappa ))), 3104 3105 3106 3107 3108 3109 3110 deriv = eval(substitute(expression({ 3111 M1 <- 2 3112 Mdiv2 <- ncol(eta) / M1 # extra$Mdiv2 3113 ymat <- matrix(y, n, Mdiv2) 3114 vTF <- c(TRUE, FALSE) 3115 locat <- eta2theta(eta[, vTF, drop = FALSE], .llocat , earg = .elocat ) 3116 Scale <- eta2theta(eta[, !vTF, drop = FALSE], .lscale , earg = .escale ) 3117 kappamat <- matrix(extra$kappa, n, Mdiv2, byrow = TRUE) 3118 zedd <- abs(ymat - locat) / Scale 3119 dl.dlocat <- sqrt(2) * ifelse(ymat >= locat, kappamat, 1/kappamat) * 3120 sign(ymat - locat) / Scale 3121 dl.dscale <- sqrt(2) * ifelse(ymat >= locat, kappamat, 1/kappamat) * 3122 zedd / Scale - 1 / Scale 3123 dlocat.deta <- dtheta.deta(locat, .llocat , earg = .elocat ) 3124 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 3125 3126 ans <- c(w) * cbind(dl.dlocat * dlocat.deta, 3127 dl.dscale * dscale.deta) 3128 ans[, interleave.VGAM(ncol(ans), M1 = M1)] 3129 }), list( .escale = escale, .lscale = lscale, 3130 .elocat = elocat, .llocat = llocat, 3131 .kappa = kappa ))), 3132 weight = eval(substitute(expression({ 3133 wz <- matrix(NA_real_, n, M) 3134 d2l.dlocat2 <- 2 / Scale^2 3135 d2l.dscale2 <- 1 / Scale^2 3136 wz[, vTF] <- d2l.dlocat2 * dlocat.deta^2 3137 wz[, !vTF] <- d2l.dscale2 * dscale.deta^2 3138 c(w) * wz 3139 }), list( .escale = escale, .lscale = lscale, 3140 .elocat = elocat, .llocat = llocat )))) 3141} # End of alaplace2(). 3142 3143 3144 3145 3146 3147 3148 3149 3150 3151 3152 3153alaplace1.control <- function(maxit = 100, ...) { 3154 list(maxit = maxit) 3155} 3156 3157 3158 3159 3160 3161 3162 3163 3164 3165 alaplace1 <- 3166 function(tau = NULL, 3167 llocation = "identitylink", 3168 ilocation = NULL, 3169 kappa = sqrt(tau/(1-tau)), 3170 Scale.arg = 1, 3171 ishrinkage = 0.95, 3172 parallel.locat = TRUE ~ 0, # FALSE, 3173 digt = 4, 3174 idf.mu = 3, 3175 zero = NULL, 3176 imethod = 1) { 3177 3178 3179 3180 apply.parint.locat <- FALSE 3181 3182 3183 3184 3185 if (!is.Numeric(kappa, positive = TRUE)) 3186 stop("bad input for argument 'kappa'") 3187 if (length(tau) && 3188 max(abs(kappa - sqrt(tau/(1-tau)))) > 1.0e-6) 3189 stop("arguments 'kappa' and 'tau' do not match") 3190 3191 if (!is.Numeric(imethod, length.arg = 1, 3192 integer.valued = TRUE, positive = TRUE) || 3193 imethod > 4) 3194 stop("argument 'imethod' must be 1, 2 or ... 4") 3195 3196 3197 llocation <- llocation 3198 3199 llocat <- as.list(substitute(llocation)) 3200 elocat <- link2list(llocat) 3201 llocat <- attr(elocat, "function.name") 3202 ilocat <- ilocation 3203 3204 3205 if (!is.Numeric(ishrinkage, length.arg = 1) || 3206 ishrinkage < 0 || 3207 ishrinkage > 1) 3208 stop("bad input for argument 'ishrinkage'") 3209 if (!is.Numeric(Scale.arg, positive = TRUE)) 3210 stop("bad input for argument 'Scale.arg'") 3211 3212 3213 3214 3215 fittedMean <- FALSE 3216 if (!is.logical(fittedMean) || length(fittedMean) != 1) 3217 stop("bad input for argument 'fittedMean'") 3218 3219 3220 3221 3222 3223 new("vglmff", 3224 blurb = c("One-parameter asymmetric Laplace distribution\n\n", 3225 "Links: ", 3226 namesof("location", llocat, earg = elocat), 3227 "\n", "\n", 3228 "Mean: location + scale * (1/kappa - kappa) / ", 3229 "sqrt(2)", "\n", 3230 "Quantiles: location", "\n", 3231 "Variance: scale^2 * (1 + kappa^4) / (2 * kappa^2)"), 3232 3233 3234 3235 3236 constraints = eval(substitute(expression({ 3237 3238 onemat <- matrix(1, M, 1) 3239 constraints.orig <- constraints 3240 3241 3242 cm1.locat <- diag(M) 3243 cmk.locat <- onemat 3244 con.locat <- cm.VGAM(cmk.locat, 3245 x = x, bool = .parallel.locat , 3246 constraints = constraints.orig, 3247 apply.int = .apply.parint.locat , 3248 cm.default = cm1.locat, 3249 cm.intercept.default = cm1.locat) 3250 3251 3252 constraints <- con.locat 3253 3254 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 3255 predictors.names = predictors.names, 3256 M1 = 1) 3257 }), list( .parallel.locat = parallel.locat, 3258 .zero = zero, 3259 .apply.parint.locat = apply.parint.locat ))), 3260 3261 3262 3263 infos = eval(substitute(function(...) { 3264 list(M1 = 1, 3265 Q1 = 1, 3266 summary.pvalues = FALSE, 3267 tau = .tau , 3268 multipleResponses = FALSE, 3269 parameters.names = c("location"), 3270 kappa = .kappa) 3271 }, list( .kappa = kappa, 3272 .tau = tau ))), 3273 initialize = eval(substitute(expression({ 3274 extra$M1 <- M1 <- 1 3275 3276 3277 temp5 <- 3278 w.y.check(w = w, y = y, 3279 ncol.w.max = if (length( .kappa ) > 1) 1 else Inf, 3280 ncol.y.max = if (length( .kappa ) > 1) 1 else Inf, 3281 out.wy = TRUE, 3282 maximize = TRUE) 3283 w <- temp5$w 3284 y <- temp5$y 3285 3286 3287 3288 extra$ncoly <- ncoly <- ncol(y) 3289 if ((ncoly > 1) && (length( .kappa ) > 1 || 3290 length( .Scale.arg ) > 1)) 3291 stop("response must be a vector if 'kappa' or 'Scale.arg' ", 3292 "has a length greater than one") 3293 3294 extra$kappa <- .kappa 3295 extra$tau <- extra$kappa^2 / (1 + extra$kappa^2) 3296 3297 3298 extra$M <- M <- max(length( .Scale.arg ), 3299 ncoly, 3300 length( .kappa )) # Recycle 3301 extra$Scale <- rep_len( .Scale.arg , M) 3302 extra$kappa <- rep_len( .kappa , M) 3303 extra$tau <- extra$kappa^2 / (1 + extra$kappa^2) 3304 extra$n <- n 3305 3306 3307 3308 3309 extra$tau.names <- tau.names <- 3310 paste("(tau = ", round(extra$tau, digits = .digt), ")", sep = "") 3311 extra$Y.names <- Y.names <- if (ncoly > 1) dimnames(y)[[2]] else "y" 3312 if (is.null(Y.names) || any(Y.names == "")) 3313 extra$Y.names <- Y.names <- paste("y", 1:ncoly, sep = "") 3314 extra$y.names <- y.names <- 3315 if (ncoly > 1) paste(Y.names, tau.names, sep = "") else tau.names 3316 3317 extra$individual <- FALSE 3318 3319 mynames1 <- param.names("location", M, skip1 = TRUE) 3320 predictors.names <- 3321 c(namesof(mynames1, .llocat , earg = .elocat , tag = FALSE)) 3322 3323 3324 locat.init <- matrix(0, n, M) 3325 if (!length(etastart)) { 3326 3327 for (jay in 1:M) { 3328 y.use <- if (ncoly > 1) y[, jay] else y 3329 if ( .imethod == 1) { 3330 locat.init[, jay] <- weighted.mean(y.use, w[, min(jay, ncol(w))]) 3331 } else if ( .imethod == 2) { 3332 locat.init[, jay] <- median(y.use) 3333 } else if ( .imethod == 3) { 3334 Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 3335 y = y.use, w = w, df = .idf.mu ) 3336 locat.init[, jay] <- c(predict(Fit5, 3337 x = x[, min(ncol(x), 2)])$y) 3338 } else { 3339 use.this <- weighted.mean(y.use, w[, min(jay, ncol(w))]) 3340 locat.init[, jay] <- (1- .ishrinkage ) * y.use + 3341 .ishrinkage * use.this 3342 } 3343 3344 3345 if (length( .ilocat )) { 3346 locat.init <- matrix( .ilocat , n, M, byrow = TRUE) 3347 } 3348 3349 if ( .llocat == "loglink") locat.init <- abs(locat.init) 3350 etastart <- 3351 cbind(theta2eta(locat.init, .llocat , earg = .elocat )) 3352 } 3353 } 3354 }), list( .imethod = imethod, 3355 .idf.mu = idf.mu, 3356 .ishrinkage = ishrinkage, .digt = digt, 3357 .elocat = elocat, .Scale.arg = Scale.arg, 3358 .llocat = llocat, .kappa = kappa, 3359 .ilocat = ilocat ))), 3360 linkinv = eval(substitute(function(eta, extra = NULL) { 3361 if ( .fittedMean ) { 3362 kappamat <- matrix(extra$kappa, extra$n, extra$M, byrow = TRUE) 3363 locat <- eta2theta(eta, .llocat , earg = .elocat ) 3364 Scale <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 3365 locat + Scale * (1/kappamat - kappamat) 3366 } else { 3367 locat <- eta2theta(eta, .llocat , earg = .elocat ) 3368 if (length(locat) > extra$n) 3369 dimnames(locat) <- list(dimnames(eta)[[1]], extra$y.names) 3370 locat 3371 } 3372 }, list( .elocat = elocat, .llocat = llocat, 3373 .fittedMean = fittedMean, .Scale.arg = Scale.arg, 3374 .kappa = kappa ))), 3375 last = eval(substitute(expression({ 3376 M1 <- extra$M1 3377 misc$M1 <- M1 3378 misc$multipleResponses <- TRUE 3379 3380 3381 3382 misc$link <- setNames(rep_len( .llocat , M), mynames1) 3383 3384 3385 3386 3387 3388 3389 3390 misc$earg <- vector("list", M) 3391 names(misc$earg) <- names(misc$link) 3392 for (ii in 1:M) { 3393 misc$earg[[ii]] <- .elocat 3394 } 3395 3396 3397 misc$expected <- TRUE 3398 extra$kappa <- misc$kappa <- .kappa 3399 extra$tau <- misc$tau <- misc$kappa^2 / (1 + misc$kappa^2) 3400 misc$true.mu <- .fittedMean # @fitted is not a true mu? 3401 3402 extra$percentile <- numeric(M) 3403 locat <- as.matrix(locat) 3404 for (ii in 1:M) { 3405 y.use <- if (ncoly > 1) y[, ii] else y 3406 extra$percentile[ii] <- 3407 100 * weighted.mean(y.use <= locat[, ii], w[, min(ii, ncol(w))]) 3408 } 3409 names(extra$percentile) <- y.names 3410 3411 extra$Scale.arg <- .Scale.arg 3412 }), list( .elocat = elocat, 3413 .llocat = llocat, 3414 .Scale.arg = Scale.arg, .fittedMean = fittedMean, 3415 .kappa = kappa ))), 3416 3417 loglikelihood = eval(substitute( 3418 function(mu, y, w, residuals = FALSE, eta, 3419 extra = NULL, 3420 summation = TRUE) { 3421 3422 ymat <- matrix(y, extra$n, extra$M) 3423 locat <- eta2theta(eta, .llocat , earg = .elocat ) 3424 kappamat <- matrix(extra$kappa, extra$n, extra$M, byrow = TRUE) 3425 Scale <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 3426 3427 if (residuals) { 3428 stop("loglikelihood residuals not implemented yet") 3429 } else { 3430 ll.elts <- c(w) * dalap(x = c(ymat), locat = c(locat), 3431 scale = c(Scale), kappa = c(kappamat), 3432 log = TRUE) 3433 if (summation) { 3434 sum(ll.elts) 3435 } else { 3436 ll.elts 3437 } 3438 } 3439 }, list( .elocat = elocat, 3440 .llocat = llocat, 3441 .Scale.arg = Scale.arg, .kappa = kappa ))), 3442 vfamily = c("alaplace1"), 3443 validparams = eval(substitute(function(eta, y, extra = NULL) { 3444 locat <- eta2theta(eta, .llocat , earg = .elocat ) 3445 okay1 <- all(is.finite(locat)) 3446 okay1 3447 }, list( .elocat = elocat, 3448 .llocat = llocat, 3449 .Scale.arg = Scale.arg, .kappa = kappa ))), 3450 3451 3452 simslot = eval(substitute( 3453 function(object, nsim) { 3454 3455 pwts <- if (length(pwts <- object@prior.weights) > 0) 3456 pwts else weights(object, type = "prior") 3457 if (any(pwts != 1)) 3458 warning("ignoring prior weights") 3459 eta <- predict(object) 3460 extra <- object@extra 3461 locat <- eta2theta(eta, .llocat , .elocat ) 3462 Scale <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 3463 kappamat <- matrix(extra$kappa, extra$n, extra$M, byrow = TRUE) 3464 ralap(nsim * length(Scale), location = c(locat), 3465 scale = c(Scale), kappa = c(kappamat)) 3466 }, list( .elocat = elocat, .llocat = llocat, 3467 .Scale.arg = Scale.arg, .kappa = kappa ))), 3468 3469 3470 3471 deriv = eval(substitute(expression({ 3472 ymat <- matrix(y, n, M) 3473 Scale <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 3474 3475 locat <- eta2theta(eta, .llocat , earg = .elocat ) 3476 3477 kappamat <- matrix(extra$kappa, n, M, byrow = TRUE) 3478 zedd <- abs(ymat-locat) / Scale 3479 3480 dl.dlocat <- ifelse(ymat >= locat, kappamat, 1/kappamat) * 3481 sqrt(2) * sign(ymat - locat) / Scale 3482 dlocat.deta <- dtheta.deta(locat, .llocat , earg = .elocat ) 3483 3484 c(w) * cbind(dl.dlocat * dlocat.deta) 3485 }), list( .Scale.arg = Scale.arg, .elocat = elocat, 3486 .llocat = llocat, .kappa = kappa ))), 3487 3488 weight = eval(substitute(expression({ 3489 d2l.dlocat2 <- 2 / Scale^2 3490 wz <- cbind(d2l.dlocat2 * dlocat.deta^2) 3491 3492 c(w) * wz 3493 }), list( .Scale.arg = Scale.arg, 3494 .elocat = elocat, .llocat = llocat )))) 3495} 3496 3497 3498 3499 3500 3501 3502 3503 3504 3505alaplace3.control <- function(maxit = 100, ...) { 3506 list(maxit = maxit) 3507} 3508 3509 3510 3511 3512 alaplace3 <- 3513 function(llocation = "identitylink", lscale = "loglink", lkappa = "loglink", 3514 ilocation = NULL, iscale = NULL, ikappa = 1.0, 3515 imethod = 1, zero = c("scale", "kappa")) { 3516 3517 llocat <- as.list(substitute(llocation)) 3518 elocat <- link2list(llocat) 3519 llocat <- attr(elocat, "function.name") 3520 ilocat <- ilocation 3521 3522 lscale <- as.list(substitute(lscale)) 3523 escale <- link2list(lscale) 3524 lscale <- attr(escale, "function.name") 3525 3526 lkappa <- as.list(substitute(lkappa)) 3527 ekappa <- link2list(lkappa) 3528 lkappa <- attr(ekappa, "function.name") 3529 3530 3531 if (!is.Numeric(imethod, length.arg = 1, 3532 integer.valued = TRUE, positive = TRUE) || 3533 imethod > 2) 3534 stop("argument 'imethod' must be 1 or 2") 3535 if (length(iscale) && 3536 !is.Numeric(iscale, positive = TRUE)) 3537 stop("bad input for argument 'iscale'") 3538 3539 3540 new("vglmff", 3541 blurb = c("Three-parameter asymmetric Laplace distribution\n\n", 3542 "Links: ", 3543 namesof("location", llocat, earg = elocat), ", ", 3544 namesof("scale", lscale, earg = escale), ", ", 3545 namesof("kappa", lkappa, earg = ekappa), 3546 "\n", "\n", 3547 "Mean: location + scale * (1/kappa - kappa) / sqrt(2)", 3548 "\n", 3549 "Variance: Scale^2 * (1 + kappa^4) / (2 * kappa^2)"), 3550 constraints = eval(substitute(expression({ 3551 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 3552 predictors.names = predictors.names, 3553 M1 = 3) 3554 }), list( .zero = zero ))), 3555 3556 infos = eval(substitute(function(...) { 3557 list(M1 = 3, 3558 Q1 = 1, 3559 multipleResponses = FALSE, 3560 parameters.names = c("location", "scale", "kappa"), 3561 summary.pvalues = FALSE, 3562 zero = .zero ) 3563 }, list( .zero = zero ))), 3564 3565 initialize = eval(substitute(expression({ 3566 3567 w.y.check(w = w, y = y, 3568 ncol.w.max = 1, 3569 ncol.y.max = 1) 3570 3571 3572 3573 predictors.names <- 3574 c(namesof("location", .llocat , earg = .elocat, tag = FALSE), 3575 namesof("scale", .lscale , earg = .escale, tag = FALSE), 3576 namesof("kappa", .lkappa , earg = .ekappa, tag = FALSE)) 3577 3578 if (!length(etastart)) { 3579 kappa.init <- if (length( .ikappa )) 3580 rep_len( .ikappa , n) else 3581 rep_len( 1.0 , n) 3582 if ( .imethod == 1) { 3583 locat.init <- median(y) 3584 scale.init <- sqrt(var(y) / 2) 3585 } else { 3586 locat.init <- y 3587 scale.init <- sqrt(sum(c(w)*abs(y-median(y ))) / (sum(w) *2)) 3588 } 3589 locat.init <- if (length( .ilocat )) 3590 rep_len( .ilocat , n) else 3591 rep_len(locat.init, n) 3592 scale.init <- if (length( .iscale )) 3593 rep_len( .iscale , n) else 3594 rep_len(scale.init, n) 3595 etastart <- 3596 cbind(theta2eta(locat.init, .llocat , earg = .elocat ), 3597 theta2eta(scale.init, .lscale , earg = .escale ), 3598 theta2eta(kappa.init, .lkappa, earg = .ekappa)) 3599 } 3600 }), list( .imethod = imethod, 3601 .elocat = elocat, .escale = escale, .ekappa = ekappa, 3602 .llocat = llocat, .lscale = lscale, .lkappa = lkappa, 3603 .ilocat = ilocat, .iscale = iscale, .ikappa = ikappa ))), 3604 linkinv = eval(substitute(function(eta, extra = NULL) { 3605 locat <- eta2theta(eta[, 1], .llocat , earg = .elocat ) 3606 Scale <- eta2theta(eta[, 2], .lscale , earg = .escale ) 3607 kappa <- eta2theta(eta[, 3], .lkappa, earg = .ekappa) 3608 locat + Scale * (1/kappa - kappa) / sqrt(2) 3609 }, list( .elocat = elocat, .llocat = llocat, 3610 .escale = escale, .lscale = lscale, 3611 .ekappa = ekappa, .lkappa = lkappa ))), 3612 last = eval(substitute(expression({ 3613 misc$link <- c(location = .llocat , 3614 scale = .lscale , 3615 kappa = .lkappa ) 3616 3617 misc$earg <- list(location = .elocat, 3618 scale = .escale, 3619 kappa = .ekappa ) 3620 3621 misc$expected = TRUE 3622 }), list( .elocat = elocat, .llocat = llocat, 3623 .escale = escale, .lscale = lscale, 3624 .ekappa = ekappa, .lkappa = lkappa ))), 3625 loglikelihood = eval(substitute( 3626 function(mu, y, w, residuals = FALSE, eta, 3627 extra = NULL, 3628 summation = TRUE) { 3629 locat <- eta2theta(eta[, 1], .llocat , earg = .elocat ) 3630 Scale <- eta2theta(eta[, 2], .lscale , earg = .escale ) 3631 kappa <- eta2theta(eta[, 3], .lkappa , earg = .ekappa ) # a matrix 3632 if (residuals) { 3633 stop("loglikelihood residuals not implemented yet") 3634 } else { 3635 ll.elts <- c(w) * dalap(x = y, locat = locat, 3636 scale = Scale, kappa = kappa, log = TRUE) 3637 if (summation) { 3638 sum(ll.elts) 3639 } else { 3640 ll.elts 3641 } 3642 } 3643 }, list( .elocat = elocat, .llocat = llocat, 3644 .escale = escale, .lscale = lscale, 3645 .ekappa = ekappa, .lkappa = lkappa ))), 3646 vfamily = c("alaplace3"), 3647 validparams = eval(substitute(function(eta, y, extra = NULL) { 3648 locat <- eta2theta(eta[, 1], .llocat , earg = .elocat ) 3649 Scale <- eta2theta(eta[, 2], .lscale , earg = .escale ) 3650 kappa <- eta2theta(eta[, 3], .lkappa , earg = .ekappa ) 3651 okay1 <- all(is.finite(locat)) && 3652 all(is.finite(Scale)) && all(0 < Scale) && 3653 all(is.finite(kappa)) && all(0 < kappa) 3654 okay1 3655 }, list( .elocat = elocat, .llocat = llocat, 3656 .escale = escale, .lscale = lscale, 3657 .ekappa = ekappa, .lkappa = lkappa ))), 3658 deriv = eval(substitute(expression({ 3659 locat <- eta2theta(eta[, 1], .llocat , earg = .elocat ) 3660 Scale <- eta2theta(eta[, 2], .lscale , earg = .escale ) 3661 kappa <- eta2theta(eta[, 3], .lkappa , earg = .ekappa ) 3662 3663 zedd <- abs(y - locat) / Scale 3664 dl.dlocat <- sqrt(2) * ifelse(y >= locat, kappa, 1/kappa) * 3665 sign(y-locat) / Scale 3666 dl.dscale <- sqrt(2) * ifelse(y >= locat, kappa, 1/kappa) * 3667 zedd / Scale - 1 / Scale 3668 dl.dkappa <- 1 / kappa - 2 * kappa / (1+kappa^2) - 3669 (sqrt(2) / Scale) * 3670 ifelse(y > locat, 1, -1/kappa^2) * abs(y-locat) 3671 3672 dlocat.deta <- dtheta.deta(locat, .llocat , earg = .elocat ) 3673 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 3674 dkappa.deta <- dtheta.deta(kappa, .lkappa, earg = .ekappa) 3675 3676 c(w) * cbind(dl.dlocat * dlocat.deta, 3677 dl.dscale * dscale.deta, 3678 dl.dkappa * dkappa.deta) 3679 }), list( .escale = escale, .lscale = lscale, 3680 .elocat = elocat, .llocat = llocat, 3681 .ekappa = ekappa, .lkappa = lkappa ))), 3682 weight = eval(substitute(expression({ 3683 d2l.dlocat2 <- 2 / Scale^2 3684 d2l.dscale2 <- 1 / Scale^2 3685 d2l.dkappa2 <- 1 / kappa^2 + 4 / (1+kappa^2)^2 3686 d2l.dkappadloc <- -sqrt(8) / ((1+kappa^2) * Scale) 3687 d2l.dkappadscale <- -(1-kappa^2) / ((1+kappa^2) * kappa * Scale) 3688 wz <- matrix(0, nrow = n, dimm(M)) 3689 wz[,iam(1, 1, M)] <- d2l.dlocat2 * dlocat.deta^2 3690 wz[,iam(2, 2, M)] <- d2l.dscale2 * dscale.deta^2 3691 wz[,iam(3, 3, M)] <- d2l.dkappa2 * dkappa.deta^2 3692 wz[,iam(1, 3, M)] <- d2l.dkappadloc * dkappa.deta * dlocat.deta 3693 wz[,iam(2, 3, M)] <- d2l.dkappadscale * dkappa.deta * dscale.deta 3694 c(w) * wz 3695 }), list( .escale = escale, .lscale = lscale, 3696 .elocat = elocat, .llocat = llocat )))) 3697} 3698 3699 3700 3701 3702 3703 3704 3705 3706 3707dlaplace <- function(x, location = 0, scale = 1, log = FALSE) { 3708 if (!is.logical(log.arg <- log) || length(log) != 1) 3709 stop("bad input for argument 'log'") 3710 rm(log) 3711 3712 3713 logdensity <- (-abs(x-location)/scale) - log(2*scale) 3714 if (log.arg) logdensity else exp(logdensity) 3715} 3716 3717 3718 3719plaplace <- function(q, location = 0, scale = 1, 3720 lower.tail = TRUE, log.p =FALSE) { 3721 zedd <- (q - location) / scale 3722 3723 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 3724 stop("bad input for argument 'lower.tail'") 3725 if (!is.logical(log.p) || length(log.p) != 1) 3726 stop("bad input for argument 'log.p'") 3727 3728 L <- max(length(q), length(location), length(scale)) 3729 if (length(q) != L) q <- rep_len(q, L) 3730 if (length(location) != L) location <- rep_len(location, L) 3731 if (length(scale) != L) scale <- rep_len(scale, L) 3732 3733 3734 if (lower.tail) { 3735 if (log.p) { 3736 ans <- ifelse(q < location, log(0.5) + zedd, 3737 log1p(- 0.5 * exp(-zedd))) 3738 } else { 3739 ans <- ifelse(q < location, 0.5 * exp(zedd), 1 - 0.5 * exp(-zedd)) 3740 } 3741 } else { 3742 if (log.p) { 3743 ans <- ifelse(q < location, log1p(- 0.5 * exp(zedd)), 3744 log(0.5) - zedd) 3745 } else { 3746 ans <- ifelse(q < location, 1 - 0.5 * exp(zedd), 0.5 * exp(-zedd)) 3747 } 3748 } 3749 ans[scale <= 0] <- NaN 3750 ans 3751} 3752 3753 3754 3755qlaplace <- function(p, location = 0, scale = 1, 3756 lower.tail = TRUE, log.p = FALSE) { 3757 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 3758 stop("bad input for argument 'lower.tail'") 3759 3760 if (!is.logical(log.p) || length(log.p) != 1) 3761 stop("bad input for argument 'log.p'") 3762 3763 3764 L <- max(length(p), length(location), length(scale)) 3765 if (length(p) != L) p <- rep_len(p, L) 3766 if (length(location) != L) location <- rep_len(location, L) 3767 if (length(scale) != L) scale <- rep_len(scale, L) 3768 3769 3770 if (lower.tail) { 3771 if (log.p) { 3772 ln.p <- p 3773 ans <- location - sign(exp(ln.p)-0.5) * scale * 3774 log(2 * ifelse(exp(ln.p) < 0.5, exp(ln.p), -expm1(ln.p))) 3775 } else { 3776 ans <- location - sign(p-0.5) * scale * 3777 log(2 * ifelse(p < 0.5, p, 1-p)) 3778 } 3779 } else { 3780 if (log.p) { 3781 ln.p <- p 3782 ans <- location - sign(0.5 - exp(ln.p)) * scale * 3783 log(2 * ifelse(-expm1(ln.p) < 0.5, -expm1(ln.p), exp(ln.p))) 3784 # ans[ln.p > 0] <- NaN 3785 } else { 3786 ans <- location - sign(0.5 - p) * scale * 3787 log(2 * ifelse(p > 0.5, 1 - p, p)) 3788 } 3789 } 3790 3791 ans[scale <= 0] <- NaN 3792 ans 3793} 3794 3795 3796 3797rlaplace <- function(n, location = 0, scale = 1) { 3798 3799 use.n <- if ((length.n <- length(n)) > 1) length.n else 3800 if (!is.Numeric(n, integer.valued = TRUE, 3801 length.arg = 1, positive = TRUE)) 3802 stop("bad input for argument 'n'") else n 3803 3804 if (!is.Numeric(scale, positive = TRUE)) 3805 stop("'scale' must be positive") 3806 3807 location <- rep_len(location, use.n) 3808 scale <- rep_len(scale, use.n) 3809 rrrr <- runif(use.n) 3810 3811 3812 3813 location - sign(rrrr - 0.5) * scale * 3814 (log(2) + ifelse(rrrr < 0.5, log(rrrr), log1p(-rrrr))) 3815} 3816 3817 3818 3819 3820 3821 laplace <- function(llocation = "identitylink", lscale = "loglink", 3822 ilocation = NULL, iscale = NULL, 3823 imethod = 1, 3824 zero = "scale") { 3825 3826 llocat <- as.list(substitute(llocation)) 3827 elocat <- link2list(llocat) 3828 llocat <- attr(elocat, "function.name") 3829 ilocat <- ilocation 3830 3831 lscale <- as.list(substitute(lscale)) 3832 escale <- link2list(lscale) 3833 lscale <- attr(escale, "function.name") 3834 3835 3836 3837 if (!is.Numeric(imethod, length.arg = 1, 3838 integer.valued = TRUE, positive = TRUE) || 3839 imethod > 3) 3840 stop("argument 'imethod' must be 1 or 2 or 3") 3841 3842 3843 if (length(iscale) && 3844 !is.Numeric(iscale, positive = TRUE)) 3845 stop("bad input for argument 'iscale'") 3846 3847 3848 new("vglmff", 3849 blurb = c("Two-parameter Laplace distribution\n\n", 3850 "Links: ", 3851 namesof("location", llocat, earg = elocat), ", ", 3852 namesof("scale", lscale, earg = escale), 3853 "\n", "\n", 3854 "Mean: location", "\n", 3855 "Variance: 2*scale^2"), 3856 constraints = eval(substitute(expression({ 3857 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 3858 predictors.names = predictors.names, 3859 M1 = 2) 3860 }), list( .zero = zero ))), 3861 3862 infos = eval(substitute(function(...) { 3863 list(M1 = 2, 3864 Q1 = 1, 3865 multipleResponses = FALSE, 3866 parameters.names = c("location", "scale"), 3867 summary.pvalues = FALSE, 3868 zero = .zero ) 3869 }, list( .zero = zero ))), 3870 3871 initialize = eval(substitute(expression({ 3872 3873 w.y.check(w = w, y = y, 3874 ncol.w.max = 1, 3875 ncol.y.max = 1) 3876 3877 3878 3879 3880 predictors.names <- 3881 c(namesof("location", .llocat , earg = .elocat, tag = FALSE), 3882 namesof("scale", .lscale , earg = .escale, tag = FALSE)) 3883 3884 3885 if (!length(etastart)) { 3886 if ( .imethod == 1) { 3887 locat.init <- median(y) 3888 scale.init <- sqrt(var(y) / 2) 3889 } else if ( .imethod == 2) { 3890 locat.init <- weighted.mean(y, w) 3891 scale.init <- sqrt(var(y) / 2) 3892 } else { 3893 locat.init <- median(y) 3894 scale.init <- sqrt(sum(c(w)*abs(y-median(y ))) / (sum(w) *2)) 3895 } 3896 locat.init <- if (length( .ilocat )) 3897 rep_len( .ilocat , n) else 3898 rep_len(locat.init, n) 3899 scale.init <- if (length( .iscale )) 3900 rep_len( .iscale , n) else 3901 rep_len(scale.init, n) 3902 etastart <- 3903 cbind(theta2eta(locat.init, .llocat , earg = .elocat ), 3904 theta2eta(scale.init, .lscale , earg = .escale )) 3905 } 3906 }), list( .imethod = imethod, 3907 .elocat = elocat, .escale = escale, 3908 .llocat = llocat, .lscale = lscale, 3909 .ilocat = ilocat, .iscale = iscale ))), 3910 linkinv = eval(substitute(function(eta, extra = NULL) { 3911 eta2theta(eta[, 1], .llocat , earg = .elocat ) 3912 }, list( .elocat = elocat, .llocat = llocat ))), 3913 last = eval(substitute(expression({ 3914 misc$link <- c(location = .llocat , scale = .lscale ) 3915 3916 misc$earg <- list(location = .elocat , scale = .escale ) 3917 3918 misc$expected <- TRUE 3919 misc$RegCondOK <- FALSE # Save this for later 3920 }), list( .escale = escale, .lscale = lscale, 3921 .elocat = elocat, .llocat = llocat ))), 3922 loglikelihood = eval(substitute( 3923 function(mu, y, w, residuals = FALSE, eta, 3924 extra = NULL, 3925 summation = TRUE) { 3926 3927 locat <- eta2theta(eta[, 1], .llocat , earg = .elocat ) 3928 Scale <- eta2theta(eta[, 2], .lscale , earg = .escale ) 3929 if (residuals) { 3930 stop("loglikelihood residuals not implemented yet") 3931 } else { 3932 ll.elts <- c(w) * dlaplace(x = y, locat = locat, 3933 scale = Scale, log = TRUE) 3934 if (summation) { 3935 sum(ll.elts) 3936 } else { 3937 ll.elts 3938 } 3939 } 3940 }, list( .escale = escale, .lscale = lscale, 3941 .elocat = elocat, .llocat = llocat ))), 3942 vfamily = c("laplace"), 3943 validparams = eval(substitute(function(eta, y, extra = NULL) { 3944 Locat <- eta2theta(eta[, 1], .llocat , earg = .elocat ) 3945 Scale <- eta2theta(eta[, 2], .lscale , earg = .escale ) 3946 okay1 <- all(is.finite(Locat)) && 3947 all(is.finite(Scale)) && all(0 < Scale) 3948 okay1 3949 }, list( .escale = escale, .lscale = lscale, 3950 .elocat = elocat, .llocat = llocat ))), 3951 deriv = eval(substitute(expression({ 3952 Locat <- eta2theta(eta[, 1], .llocat , earg = .elocat ) 3953 Scale <- eta2theta(eta[, 2], .lscale , earg = .escale ) 3954 3955 zedd <- abs(y-Locat) / Scale 3956 dl.dLocat <- sign(y - Locat) / Scale 3957 dl.dscale <- zedd / Scale - 1 / Scale 3958 3959 dLocat.deta <- dtheta.deta(Locat, .llocat , earg = .elocat ) 3960 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 3961 3962 c(w) * cbind(dl.dLocat * dLocat.deta, 3963 dl.dscale * dscale.deta) 3964 }), list( .escale = escale, .lscale = lscale, 3965 .elocat = elocat, .llocat = llocat ))), 3966 weight = eval(substitute(expression({ 3967 d2l.dLocat2 <- d2l.dscale2 <- 1 / Scale^2 3968 wz <- matrix(0, nrow = n, ncol = M) # diagonal 3969 wz[,iam(1, 1, M)] <- d2l.dLocat2 * dLocat.deta^2 3970 wz[,iam(2, 2, M)] <- d2l.dscale2 * dscale.deta^2 3971 c(w) * wz 3972 }), list( .escale = escale, .lscale = lscale, 3973 .elocat = elocat, .llocat = llocat )))) 3974} 3975 3976 3977 3978fff.control <- function(save.weights = TRUE, ...) { 3979 list(save.weights = save.weights) 3980} 3981 3982 3983 fff <- function(link = "loglink", 3984 idf1 = NULL, idf2 = NULL, nsimEIM = 100, # ncp = 0, 3985 imethod = 1, zero = NULL) { 3986 link <- as.list(substitute(link)) 3987 earg <- link2list(link) 3988 link <- attr(earg, "function.name") 3989 3990 3991 if (!is.Numeric(imethod, length.arg = 1, 3992 integer.valued = TRUE, positive = TRUE) || 3993 imethod > 2) 3994 stop("argument 'imethod' must be 1 or 2") 3995 3996 3997 if (!is.Numeric(nsimEIM, length.arg = 1, 3998 integer.valued = TRUE) || 3999 nsimEIM <= 10) 4000 stop("argument 'nsimEIM' should be an integer greater than 10") 4001 4002 ncp <- 0 4003 if (any(ncp != 0)) 4004 warning("not sure about ncp != 0 wrt dl/dtheta") 4005 4006 4007 4008 new("vglmff", 4009 blurb = c("F-distribution\n\n", 4010 "Links: ", 4011 namesof("df1", link, earg = earg), ", ", 4012 namesof("df2", link, earg = earg), 4013 "\n", "\n", 4014 "Mean: df2/(df2-2) provided df2>2 and ncp = 0", "\n", 4015 "Variance: ", 4016 "2*df2^2*(df1+df2-2)/(df1*(df2-2)^2*(df2-4)) ", 4017 "provided df2>4 and ncp = 0"), 4018 constraints = eval(substitute(expression({ 4019 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 4020 predictors.names = predictors.names, 4021 M1 = 2) 4022 }), list( .zero = zero ))), 4023 4024 infos = eval(substitute(function(...) { 4025 list(M1 = 2, 4026 Q1 = 1, 4027 multipleResponses = FALSE, 4028 parameters.names = c("df1", "df2"), 4029 zero = .zero ) 4030 }, list( .zero = zero ))), 4031 4032 initialize = eval(substitute(expression({ 4033 4034 w.y.check(w = w, y = y, 4035 ncol.w.max = 1, 4036 ncol.y.max = 1) 4037 4038 4039 4040 predictors.names <- 4041 c(namesof("df1", .link , earg = .earg , tag = FALSE), 4042 namesof("df2", .link , earg = .earg , tag = FALSE)) 4043 4044 4045 if (!length(etastart)) { 4046 if ( .imethod == 1) { 4047 df2.init <- b <- 2*mean(y) / (mean(y)-1) 4048 df1.init <- 2*b^2*(b-2)/(var(y)*(b-2)^2 * (b-4) - 2*b^2) 4049 if (df2.init < 4) df2.init <- 5 4050 if (df1.init < 2) df1.init <- 3 4051 } else { 4052 df2.init <- b <- 2*median(y) / (median(y)-1) 4053 summy <- summary(y) 4054 var.est <- summy[5] - summy[2] 4055 df1.init <- 2*b^2*(b-2)/(var.est*(b-2)^2 * (b-4) - 2*b^2) 4056 } 4057 df1.init <- if (length( .idf1 )) 4058 rep_len( .idf1 , n) else 4059 rep_len(df1.init, n) 4060 df2.init <- if (length( .idf2 )) 4061 rep_len( .idf2 , n) else 4062 rep_len(1, n) 4063 etastart <- cbind(theta2eta(df1.init, .link , earg = .earg ), 4064 theta2eta(df2.init, .link , earg = .earg )) 4065 } 4066 }), list( .imethod = imethod, .idf1 = idf1, .earg = earg, 4067 .idf2 = idf2, .link = link ))), 4068 linkinv = eval(substitute(function(eta, extra = NULL) { 4069 df2 <- eta2theta(eta[, 2], .link , earg = .earg ) 4070 ans <- df2 * NA 4071 ans[df2 > 2] <- df2[df2 > 2] / (df2[df2 > 2] - 2) 4072 ans 4073 }, list( .link = link, .earg = earg ))), 4074 last = eval(substitute(expression({ 4075 misc$link <- c(df1 = .link , df2 = .link ) 4076 4077 misc$earg <- list(df1 = .earg , df2 = .earg ) 4078 4079 misc$nsimEIM <- .nsimEIM 4080 misc$ncp <- .ncp 4081 }), list( .link = link, .earg = earg, 4082 .ncp = ncp, 4083 .nsimEIM = nsimEIM ))), 4084 loglikelihood = eval(substitute( 4085 function(mu, y, w, residuals = FALSE, eta, 4086 extra = NULL, 4087 summation = TRUE) { 4088 4089 df1 <- eta2theta(eta[, 1], .link , earg = .earg ) 4090 df2 <- eta2theta(eta[, 2], .link , earg = .earg ) 4091 if (residuals) { 4092 stop("loglikelihood residuals not implemented yet") 4093 } else { 4094 ll.elts <- c(w) * df(x = y, df1 = df1, df2 = df2, 4095 ncp = .ncp , log = TRUE) 4096 if (summation) { 4097 sum(ll.elts) 4098 } else { 4099 ll.elts 4100 } 4101 } 4102 }, list( .link = link, .earg = earg, .ncp = ncp ))), 4103 vfamily = c("fff"), 4104 validparams = eval(substitute(function(eta, y, extra = NULL) { 4105 df1 <- eta2theta(eta[, 1], .link , earg = .earg ) 4106 df2 <- eta2theta(eta[, 2], .link , earg = .earg ) 4107 okay1 <- all(is.finite(df1)) && all(0 < df1) && 4108 all(is.finite(df2)) && all(0 < df2) 4109 okay1 4110 }, list( .link = link, .earg = earg, .ncp = ncp ))), 4111 deriv = eval(substitute(expression({ 4112 df1 <- eta2theta(eta[, 1], .link , earg = .earg ) 4113 df2 <- eta2theta(eta[, 2], .link , earg = .earg ) 4114 dl.ddf1 <- 0.5*digamma(0.5*(df1+df2)) + 0.5 + 0.5*log(df1/df2) + 4115 0.5*log(y) - 0.5*digamma(0.5*df1) - 4116 0.5*(df1+df2)*(y/df2) / (1 + df1*y/df2) - 4117 0.5*log1p(df1*y/df2) 4118 dl.ddf2 <- 0.5*digamma(0.5*(df1+df2)) - 0.5*df1/df2 - 4119 0.5*digamma(0.5*df2) - 4120 0.5*(df1+df2) * (-df1*y/df2^2) / (1 + df1*y/df2) - 4121 0.5*log1p(df1*y/df2) 4122 ddf1.deta <- dtheta.deta(df1, .link , earg = .earg ) 4123 ddf2.deta <- dtheta.deta(df2, .link , earg = .earg ) 4124 dthetas.detas <- cbind(ddf1.deta, ddf2.deta) 4125 c(w) * dthetas.detas * cbind(dl.ddf1, dl.ddf2) 4126 }), list( .link = link, .earg = earg ))), 4127 weight = eval(substitute(expression({ 4128 run.varcov <- 0 4129 ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE) 4130 for (ii in 1:( .nsimEIM )) { 4131 ysim <- rf(n = n, df1=df1, df2=df2) 4132 dl.ddf1 <- 0.5*digamma(0.5*(df1+df2)) + 0.5 + 0.5*log(df1/df2) + 4133 0.5*log(ysim) - 0.5*digamma(0.5*df1) - 4134 0.5*(df1+df2)*(ysim/df2) / (1 + df1*ysim/df2) - 4135 0.5*log1p(df1*ysim/df2) 4136 dl.ddf2 <- 0.5*digamma(0.5*(df1+df2)) - 0.5*df1/df2 - 4137 0.5*digamma(0.5*df2) - 4138 0.5*(df1+df2) * (-df1*ysim/df2^2)/(1 + df1*ysim/df2) - 4139 0.5*log1p(df1*ysim/df2) 4140 rm(ysim) 4141 temp3 <- cbind(dl.ddf1, dl.ddf2) 4142 run.varcov <- ((ii-1) * run.varcov + 4143 temp3[,ind1$row.index]*temp3[,ind1$col.index]) / ii 4144 } 4145 wz <- if (intercept.only) 4146 matrix(colMeans(run.varcov), 4147 n, ncol(run.varcov), byrow = TRUE) else run.varcov 4148 4149 wz <- c(w) * wz * dthetas.detas[, ind1$row] * 4150 dthetas.detas[, ind1$col] 4151 wz 4152 }), list( .link = link, .earg = earg, .nsimEIM = nsimEIM, 4153 .ncp = ncp )))) 4154} 4155 4156 4157 4158 4159 hyperg <- function(N = NULL, D = NULL, 4160 lprob = "logitlink", 4161 iprob = NULL) { 4162 4163 inputN <- is.Numeric(N, positive = TRUE) 4164 inputD <- is.Numeric(D, positive = TRUE) 4165 if (inputD && inputN) 4166 stop("only one of 'N' and 'D' is to be inputted") 4167 if (!inputD && !inputN) 4168 stop("one of 'N' and 'D' needs to be inputted") 4169 4170 4171 lprob <- as.list(substitute(lprob)) 4172 earg <- link2list(lprob) 4173 lprob <- attr(earg, "function.name") 4174 4175 4176 4177 new("vglmff", 4178 blurb = c("Hypergeometric distribution\n\n", 4179 "Link: ", 4180 namesof("prob", lprob, earg = earg), "\n", 4181 "Mean: D/N\n"), 4182 initialize = eval(substitute(expression({ 4183 NCOL <- function (x) 4184 if (is.array(x) && length(dim(x)) > 1 || 4185 is.data.frame(x)) ncol(x) else as.integer(1) 4186 if (NCOL(y) == 1) { 4187 if (is.factor(y)) y <- y != levels(y)[1] 4188 nn <- rep_len(1, n) 4189 if (!all(y >= 0 & y <= 1)) 4190 stop("response values must be in [0, 1]") 4191 mustart <- (0.5 + w * y) / (1 + w) 4192 no.successes <- w * y 4193 if (any(abs(no.successes - round(no.successes)) > 0.001)) 4194 stop("Number of successes must be integer-valued") 4195 } else if (NCOL(y) == 2) { 4196 if (any(abs(y - round(y)) > 0.001)) 4197 stop("Count data must be integer-valued") 4198 nn <- y[, 1] + y[, 2] 4199 y <- ifelse(nn > 0, y[, 1]/nn, 0) 4200 w <- w * nn 4201 mustart <- (0.5 + nn * y) / (1 + nn) 4202 mustart[mustart >= 1] <- 0.95 4203 } else 4204 stop("Response not of the right form") 4205 4206 predictors.names <- 4207 namesof("prob", .lprob , earg = .earg , tag = FALSE) 4208 extra$Nvector <- .N 4209 extra$Dvector <- .D 4210 extra$Nunknown <- length(extra$Nvector) == 0 4211 if (!length(etastart)) { 4212 init.prob <- if (length( .iprob)) 4213 rep_len( .iprob, n) else 4214 mustart 4215 etastart <- matrix(init.prob, n, NCOL(y)) 4216 4217 } 4218 }), list( .lprob = lprob, .earg = earg, .N = N, .D = D, 4219 .iprob = iprob ))), 4220 linkinv = eval(substitute(function(eta, extra = NULL) { 4221 eta2theta(eta, .lprob , earg = .earg ) 4222 }, list( .lprob = lprob, .earg = earg ))), 4223 last = eval(substitute(expression({ 4224 misc$link <- c("prob" = .lprob) 4225 4226 misc$earg <- list("prob" = .earg ) 4227 4228 misc$Dvector <- .D 4229 misc$Nvector <- .N 4230 }), list( .N = N, .D = D, .lprob = lprob, .earg = earg ))), 4231 linkfun = eval(substitute(function(mu, extra = NULL) { 4232 theta2eta(mu, .lprob, earg = .earg ) 4233 }, list( .lprob = lprob, .earg = earg ))), 4234 loglikelihood = eval(substitute( 4235 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 4236 summation = TRUE) { 4237 N <- extra$Nvector 4238 Dvec <- extra$Dvector 4239 prob <- mu 4240 yvec <- w * y 4241 if (residuals) { 4242 stop("loglikelihood residuals not implemented yet") 4243 } else { 4244 ll.elts <- 4245 if (extra$Nunknown) { 4246 tmp12 <- Dvec * (1-prob) / prob 4247 4248 4249 (lgamma(1+tmp12) + lgamma(1+Dvec/prob-w) - 4250 lgamma(1+tmp12-w+yvec) - lgamma(1+Dvec/prob)) 4251 } else { 4252 4253 4254 (lgamma(1+N*prob) + lgamma(1+N*(1-prob)) - 4255 lgamma(1+N*prob-yvec) - 4256 lgamma(1+N*(1-prob) -w + yvec)) 4257 } 4258 if (summation) { 4259 sum(ll.elts) 4260 } else { 4261 ll.elts 4262 } 4263 } 4264 }, list( .lprob = lprob, .earg = earg ))), 4265 vfamily = c("hyperg"), 4266 validparams = eval(substitute(function(eta, y, extra = NULL) { 4267 prob <- eta2theta(eta, .lprob , earg = .earg ) 4268 okay1 <- all(is.finite(prob)) && all(0 < prob & prob < 1) 4269 okay1 4270 }, list( .lprob = lprob, .earg = earg ))), 4271 deriv = eval(substitute(expression({ 4272 prob <- mu # equivalently, eta2theta(eta, .lprob, earg = .earg ) 4273 dprob.deta <- dtheta.deta(prob, .lprob, earg = .earg ) 4274 Dvec <- extra$Dvector 4275 Nvec <- extra$Nvector 4276 yvec <- w * y 4277 if (extra$Nunknown) { 4278 tmp72 <- -Dvec / prob^2 4279 tmp12 <- Dvec * (1-prob) / prob 4280 dl.dprob <- tmp72 * (digamma(1 + tmp12) + 4281 digamma(1 + Dvec/prob -w) - 4282 digamma(1 + tmp12-w+yvec) - digamma(1 + Dvec/prob)) 4283 } else { 4284 dl.dprob <- Nvec * (digamma(1+Nvec*prob) - 4285 digamma(1+Nvec*(1-prob)) - 4286 digamma(1+Nvec*prob-yvec) + 4287 digamma(1+Nvec*(1-prob)-w+yvec)) 4288 } 4289 c(w) * dl.dprob * dprob.deta 4290 }), list( .lprob = lprob, .earg = earg ))), 4291 weight = eval(substitute(expression({ 4292 if (extra$Nunknown) { 4293 tmp722 <- tmp72^2 4294 tmp13 <- 2*Dvec / prob^3 4295 d2l.dprob2 <- tmp722 * (trigamma(1 + tmp12) + 4296 trigamma(1 + Dvec/prob - w) - 4297 trigamma(1 + tmp12 - w + yvec) - 4298 trigamma(1 + Dvec/prob)) + 4299 tmp13 * (digamma(1 + tmp12) + 4300 digamma(1 + Dvec/prob - w) - 4301 digamma(1 + tmp12 - w + yvec) - 4302 digamma(1 + Dvec/prob)) 4303 } else { 4304 d2l.dprob2 <- Nvec^2 * (trigamma(1+Nvec*prob) + 4305 trigamma(1+Nvec*(1-prob)) - 4306 trigamma(1+Nvec*prob-yvec) - 4307 trigamma(1+Nvec*(1-prob)-w+yvec)) 4308 } 4309 d2prob.deta2 <- d2theta.deta2(prob, .lprob , earg = .earg ) 4310 4311 wz <- -(dprob.deta^2) * d2l.dprob2 4312 wz <- c(w) * wz 4313 wz[wz < .Machine$double.eps] <- .Machine$double.eps 4314 wz 4315 }), list( .lprob = lprob, .earg = earg )))) 4316} 4317 4318 4319 4320 4321 4322 4323dbenini <- function(x, y0, shape, log = FALSE) { 4324 if (!is.logical(log.arg <- log) || length(log) != 1) 4325 stop("bad input for argument 'log'") 4326 rm(log) 4327 4328 4329 N <- max(length(x), length(shape), length(y0)) 4330 if (length(x) != N) x <- rep_len(x, N) 4331 if (length(shape) != N) shape <- rep_len(shape, N) 4332 if (length(y0) != N) y0 <- rep_len(y0, N) 4333 4334 logdensity <- rep_len(log(0), N) 4335 xok <- (x > y0) 4336 tempxok <- log(x[xok]/y0[xok]) 4337 logdensity[xok] <- log(2*shape[xok]) - shape[xok] * tempxok^2 + 4338 log(tempxok) - log(x[xok]) 4339 logdensity[is.infinite(x)] <- log(0) # 20141209 KaiH 4340 if (log.arg) logdensity else exp(logdensity) 4341} 4342 4343 4344 4345pbenini <- function(q, y0, shape, lower.tail = TRUE, log.p = FALSE) { 4346 if (!is.Numeric(q)) 4347 stop("bad input for argument 'q'") 4348 if (!is.Numeric(shape, positive = TRUE)) 4349 stop("bad input for argument 'shape'") 4350 if (!is.Numeric(y0, positive = TRUE)) 4351 stop("bad input for argument 'y0'") 4352 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 4353 stop("bad input for argument 'lower.tail'") 4354 if (!is.logical(log.p) || length(log.p) != 1) 4355 stop("bad input for argument 'log.p'") 4356 4357 N <- max(length(q), length(shape), length(y0)) 4358 if (length(q) != N) q <- rep_len(q, N) 4359 if (length(shape) != N) shape <- rep_len(shape, N) 4360 if (length(y0) != N) y0 <- rep_len(y0, N) 4361 4362 ans <- y0 * 0 4363 ok <- q > y0 4364 4365 4366 if (lower.tail) { 4367 if (log.p) { 4368 ans[ok] <- log(-expm1(-shape[ok] * (log(q[ok]/y0[ok]))^2)) 4369 ans[q <= y0 ] <- -Inf 4370 } else { 4371 ans[ok] <- -expm1(-shape[ok] * (log(q[ok]/y0[ok]))^2) 4372 } 4373 } else { 4374 if (log.p) { 4375 ans[ok] <- -shape[ok] * (log(q[ok]/y0[ok]))^2 4376 ans[q <= y0] <- 0 4377 } else { 4378 ans[ok] <- exp(-shape[ok] * (log(q[ok]/y0[ok]))^2) 4379 ans[q <= y0] <- 1 4380 } 4381 } 4382 4383 ans 4384} 4385 4386 4387 4388qbenini <- function(p, y0, shape, lower.tail = TRUE, log.p = FALSE) { 4389 4390 4391 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 4392 stop("bad input for argument 'lower.tail'") 4393 if (!is.logical(log.p) || length(log.p) != 1) 4394 stop("bad input for argument 'log.p'") 4395 4396 if (lower.tail) { 4397 if (log.p) { 4398 ln.p <- p 4399 ans <- y0 * exp(sqrt(-log(-expm1(ln.p)) / shape)) 4400 } else { 4401 ans <- y0 * exp(sqrt(-log1p(-p) / shape)) 4402 } 4403 } else { 4404 if (log.p) { 4405 ln.p <- p 4406 ans <- y0 * exp(sqrt(-ln.p / shape)) 4407 } else { 4408 ans <- y0 * exp(sqrt(-log(p) / shape)) 4409 } 4410 } 4411 ans[y0 <= 0] <- NaN 4412 ans 4413} 4414 4415 4416 4417rbenini <- function(n, y0, shape) { 4418 y0 * exp(sqrt(-log(runif(n)) / shape)) 4419} 4420 4421 4422 4423 4424 4425 4426 4427 benini1 <- 4428 function(y0 = stop("argument 'y0' must be specified"), 4429 lshape = "loglink", 4430 ishape = NULL, imethod = 1, zero = NULL, 4431 parallel = FALSE, 4432 type.fitted = c("percentiles", "Qlink"), 4433 percentiles = 50) { 4434 4435 type.fitted <- match.arg(type.fitted, 4436 c("percentiles", "Qlink"))[1] 4437 4438 4439 lshape <- as.list(substitute(lshape)) 4440 eshape <- link2list(lshape) 4441 lshape <- attr(eshape, "function.name") 4442 4443 4444 if (!is.Numeric(imethod, length.arg = 1, 4445 integer.valued = TRUE, positive = TRUE) || 4446 imethod > 2) 4447 stop("argument 'imethod' must be 1 or 2") 4448 4449 4450 if (!is.Numeric(y0, positive = TRUE, length.arg = 1)) 4451 stop("bad input for argument 'y0'") 4452 4453 4454 4455 4456 4457 new("vglmff", 4458 blurb = c("1-parameter Benini distribution\n\n", 4459 "Link: ", 4460 namesof("shape", lshape, earg = eshape), 4461 "\n", "\n", 4462 "Median: qbenini(p = 0.5, y0, shape)"), 4463 constraints = eval(substitute(expression({ 4464 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 4465 bool = .parallel , 4466 constraints, apply.int = FALSE) 4467 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 4468 predictors.names = predictors.names, 4469 M1 = 1) 4470 }), list( .parallel = parallel, 4471 .zero = zero ))), 4472 4473 infos = eval(substitute(function(...) { 4474 list(M1 = 1, 4475 Q1 = 1, 4476 expected = TRUE, 4477 multipleResponses = TRUE, 4478 parameters.names = c("shape"), 4479 parallel = .parallel , 4480 percentiles = .percentiles , 4481 type.fitted = .type.fitted , 4482 lshape = .lshape , 4483 eshape = .eshape , 4484 zero = .zero ) 4485 }, list( .parallel = parallel, 4486 .zero = zero, 4487 .percentiles = percentiles , 4488 .type.fitted = type.fitted, 4489 .eshape = eshape, 4490 .lshape = lshape))), 4491 4492 initialize = eval(substitute(expression({ 4493 4494 temp5 <- 4495 w.y.check(w = w, y = y, 4496 ncol.w.max = Inf, 4497 ncol.y.max = Inf, 4498 out.wy = TRUE, 4499 colsyperw = 1, 4500 maximize = TRUE) 4501 w <- temp5$w 4502 y <- temp5$y 4503 4504 4505 4506 ncoly <- ncol(y) 4507 M1 <- 1 4508 M <- M1 * ncoly 4509 extra$ncoly <- ncoly 4510 extra$type.fitted <- .type.fitted 4511 extra$colnames.y <- colnames(y) 4512 extra$percentiles <- .percentiles 4513 extra$M1 <- M1 4514 4515 4516 mynames1 <- paste("shape", if (ncoly > 1) 1:ncoly else "", sep = "") 4517 predictors.names <- 4518 namesof(mynames1, .lshape , earg = .eshape , tag = FALSE) 4519 4520 extra$y0 <- .y0 # Of unit length; 20181205; to make things easy. 4521 if (any(y <= extra$y0)) 4522 stop("some values of the response are > argument 'y0' values") 4523 4524 4525 if (!length(etastart)) { 4526 probs.y <- (1:3) / 4 4527 qofy <- quantile(rep(y, times = w), probs = probs.y) 4528 if ( .imethod == 1) { 4529 shape.init <- mean(-log1p(-probs.y) / (log(qofy))^2) 4530 } else { 4531 shape.init <- median(-log1p(-probs.y) / (log(qofy))^2) 4532 } 4533 shape.init <- matrix(if (length( .ishape )) .ishape else shape.init, 4534 n, ncoly, byrow = TRUE) 4535 etastart <- cbind(theta2eta(shape.init, .lshape , earg = .eshape )) 4536 } 4537 }), list( .imethod = imethod, 4538 .ishape = ishape, 4539 .lshape = lshape, .eshape = eshape, 4540 .percentiles = percentiles, 4541 .type.fitted = type.fitted, 4542 .y0 = y0 ))), 4543 4544 4545 4546 linkinv = eval(substitute(function(eta, extra = NULL) { 4547 type.fitted <- 4548 if (length(extra$type.fitted)) { 4549 extra$type.fitted 4550 } else { 4551 warning("cannot find 'type.fitted'. Returning the 'median'.") 4552 extra$percentiles <- 50 # Overwrite whatever was there 4553 "percentiles" 4554 } 4555 type.fitted <- match.arg(type.fitted, 4556 c("percentiles", "Qlink"))[1] 4557 4558 if (type.fitted == "Qlink") { 4559 eta2theta(eta, link = "loglink") 4560 } else { 4561 shape <- eta2theta(eta, .lshape , earg = .eshape ) 4562 4563 pcent <- extra$percentiles 4564 perc.mat <- matrix(pcent, NROW(eta), length(pcent), 4565 byrow = TRUE) / 100 4566 fv <- 4567 switch(type.fitted, 4568 "percentiles" = qbenini(perc.mat, 4569 y0 = extra$y0, 4570 shape = matrix(shape, nrow(perc.mat), ncol(perc.mat)))) 4571 if (type.fitted == "percentiles") 4572 fv <- label.cols.y(fv, colnames.y = extra$colnames.y, 4573 NOS = NCOL(eta), percentiles = pcent, 4574 one.on.one = FALSE) 4575 fv 4576 } 4577 }, list( .lshape = lshape, .eshape = eshape ))), 4578 4579 4580 4581 4582 4583 4584 4585 4586 4587 4588 last = eval(substitute(expression({ 4589 M1 <- extra$M1 4590 misc$link <- c(rep_len( .lshape , ncoly)) 4591 names(misc$link) <- mynames1 4592 4593 misc$earg <- vector("list", M) 4594 names(misc$earg) <- mynames1 4595 for (ii in 1:ncoly) { 4596 misc$earg[[ii]] <- .eshape 4597 } 4598 4599 extra$y0 <- .y0 4600 4601 }), list( .lshape = lshape, 4602 .eshape = eshape, .y0 = y0 ))), 4603 loglikelihood = eval(substitute( 4604 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 4605 summation = TRUE) { 4606 shape <- eta2theta(eta, .lshape , earg = .eshape ) 4607 y0 <- extra$y0 4608 if (residuals) { 4609 stop("loglikelihood residuals not implemented yet") 4610 } else { 4611 ll.elts <- c(w) * dbenini(x = y, y0 = y0, shape = shape, log = TRUE) 4612 if (summation) { 4613 sum(ll.elts) 4614 } else { 4615 ll.elts 4616 } 4617 } 4618 }, list( .lshape = lshape, .eshape = eshape ))), 4619 vfamily = c("benini1"), 4620 validparams = eval(substitute(function(eta, y, extra = NULL) { 4621 shape <- eta2theta(eta, .lshape , earg = .eshape ) 4622 okay1 <- all(is.finite(shape)) && all(0 < shape) 4623 okay1 4624 }, list( .lshape = lshape, .eshape = eshape ))), 4625 4626 4627 4628 4629 4630 simslot = eval(substitute( 4631 function(object, nsim) { 4632 4633 pwts <- if (length(pwts <- object@prior.weights) > 0) 4634 pwts else weights(object, type = "prior") 4635 if (any(pwts != 1)) 4636 warning("ignoring prior weights") 4637 eta <- predict(object) 4638 extra <- object@extra 4639 shape <- eta2theta(eta, .lshape , earg = .eshape ) 4640 y0 <- extra$y0 4641 rbenini(nsim * length(shape), y0 = y0, shape = shape) 4642 }, list( .lshape = lshape, .eshape = eshape ))), 4643 4644 4645 4646 4647 4648 deriv = eval(substitute(expression({ 4649 shape <- eta2theta(eta, .lshape , earg = .eshape ) 4650 4651 y0 <- extra$y0 4652 dl.dshape <- 1/shape - (log(y/y0))^2 4653 4654 dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape ) 4655 4656 c(w) * dl.dshape * dshape.deta 4657 }), list( .lshape = lshape, .eshape = eshape ))), 4658 weight = eval(substitute(expression({ 4659 ned2l.dshape2 <- 1 / shape^2 4660 wz <- ned2l.dshape2 * dshape.deta^2 4661 c(w) * wz 4662 }), list( .lshape = lshape, .eshape = eshape )))) 4663} 4664 4665 4666 4667 4668 4669 4670 dpolono <- function (x, meanlog = 0, sdlog = 1, bigx = 170, ...) { 4671 mapply(function(x, meanlog, sdlog, ...) { 4672 if (abs(x) > floor(x)) { # zero prob for -ve or non-integer 4673 0 4674 } else 4675 if (x == Inf) { # 20141215 KaiH 4676 0 4677 } else 4678 if (x > bigx) { 4679 z <- (log(x) - meanlog) / sdlog 4680 (1 + (z^2 + log(x) - meanlog - 1) / (2 * x * sdlog^2)) * 4681 exp(-0.5 * z^2) / (sqrt(2 * pi) * sdlog * x) 4682 } else 4683 integrate( function(t) exp(t * x - exp(t) - 4684 0.5 * ((t - meanlog) / sdlog)^2), 4685 lower = -Inf, upper = Inf, ...)$value / (sqrt(2 * pi) * 4686 sdlog * exp(lgamma(x + 1.0))) 4687 }, x, meanlog, sdlog, ...) 4688} 4689 4690 4691 4692 4693ppolono <- function(q, meanlog = 0, sdlog = 1, 4694 isOne = 1 - sqrt( .Machine$double.eps ), ...) { 4695 4696 4697 .cumprob <- rep_len(0, length(q)) 4698 .cumprob[q == Inf] <- 1 # special case 4699 4700 4701 q <- floor(q) 4702 ii <- -1 4703 while (any(xActive <- ((.cumprob < isOne) & (q > ii)))) 4704 .cumprob[xActive] <- .cumprob[xActive] + 4705 dpolono(ii <- (ii+1), meanlog, sdlog, ...) 4706 .cumprob 4707} 4708 4709 4710 4711 4712 4713 4714 4715 4716 4717rpolono <- function(n, meanlog = 0, sdlog = 1) { 4718 lambda <- rlnorm(n = n, meanlog = meanlog, sdlog = sdlog) 4719 rpois(n = n, lambda = lambda) 4720} 4721 4722 4723 4724 4725 4726 4727 4728 4729 4730 4731 4732 4733 4734 4735dtriangle <- function(x, theta, lower = 0, upper = 1, log = FALSE) { 4736 if (!is.logical(log.arg <- log) || length(log) != 1) 4737 stop("bad input for argument 'log'") 4738 rm(log) 4739 4740 4741 N <- max(length(x), length(theta), length(lower), length(upper)) 4742 if (length(x) != N) x <- rep_len(x, N) 4743 if (length(theta) != N) theta <- rep_len(theta, N) 4744 if (length(lower) != N) lower <- rep_len(lower, N) 4745 if (length(upper) != N) upper <- rep_len(upper, N) 4746 4747 denom1 <- ((upper-lower)*(theta-lower)) 4748 denom2 <- ((upper-lower)*(upper-theta)) 4749 logdensity <- rep_len(log(0), N) 4750 xok.neg <- (lower < x) & (x <= theta) 4751 xok.pos <- (theta <= x) & (x < upper) 4752 logdensity[xok.neg] = 4753 log(2 * (x[xok.neg] - lower[xok.neg]) / denom1[xok.neg]) 4754 logdensity[xok.pos] = 4755 log(2 * (upper[xok.pos] - x[xok.pos]) / denom2[xok.pos]) 4756 logdensity[lower >= upper] <- NaN 4757 logdensity[lower > theta] <- NaN 4758 logdensity[upper < theta] <- NaN 4759 if (log.arg) logdensity else exp(logdensity) 4760} 4761 4762 4763rtriangle <- function(n, theta, lower = 0, upper = 1) { 4764 4765 4766 use.n <- if ((length.n <- length(n)) > 1) length.n else 4767 if (!is.Numeric(n, integer.valued = TRUE, 4768 length.arg = 1, positive = TRUE)) 4769 stop("bad input for argument 'n'") else n 4770 4771 4772 if (!is.Numeric(theta)) 4773 stop("bad input for argument 'theta'") 4774 if (!is.Numeric(lower)) 4775 stop("bad input for argument 'lower'") 4776 if (!is.Numeric(upper)) 4777 stop("bad input for argument 'upper'") 4778 if (!all(lower < theta & theta < upper)) 4779 stop("lower < theta < upper values are required") 4780 4781 N <- use.n 4782 lower <- rep_len(lower, N) 4783 upper <- rep_len(upper, N) 4784 theta <- rep_len(theta, N) 4785 t1 <- sqrt(runif(n)) 4786 t2 <- sqrt(runif(n)) 4787 ifelse(runif(n) < (theta - lower) / (upper - lower), 4788 lower + (theta - lower) * t1, 4789 upper - (upper - theta) * t2) 4790} 4791 4792 4793 4794qtriangle <- function(p, theta, lower = 0, upper = 1, 4795 lower.tail = TRUE, log.p = FALSE) { 4796 4797 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 4798 stop("bad input for argument 'lower.tail'") 4799 4800 if (!is.logical(log.p) || length(log.p) != 1) 4801 stop("bad input for argument 'log.p'") 4802 4803 N <- max(length(p), length(theta), length(lower), length(upper)) 4804 if (length(p) != N) p <- rep_len(p, N) 4805 if (length(theta) != N) theta <- rep_len(theta, N) 4806 if (length(lower) != N) lower <- rep_len(lower, N) 4807 if (length(upper) != N) upper <- rep_len(upper, N) 4808 4809 ans <- NA_real_ * p 4810 if (lower.tail) { 4811 if (log.p) { 4812 Neg <- (exp(ln.p) <= (theta - lower) / (upper - lower)) 4813 temp1 <- exp(ln.p) * (upper - lower) * (theta - lower) 4814 Pos <- (exp(ln.p) >= (theta - lower) / (upper - lower)) 4815 pstar <- (exp(ln.p) - (theta - lower) / (upper - lower)) / 4816 ((upper - theta) / (upper - lower)) 4817 } else { 4818 Neg <- (p <= (theta - lower) / (upper - lower)) 4819 temp1 <- p * (upper - lower) * (theta - lower) 4820 Pos <- (p >= (theta - lower) / (upper - lower)) 4821 pstar <- (p - (theta - lower) / (upper - lower)) / 4822 ((upper - theta) / (upper - lower)) 4823 } 4824 } else { 4825 if (log.p) { 4826 ln.p <- p 4827 Neg <- (exp(ln.p) >= (upper- theta) / (upper - lower)) 4828 temp1 <- -expm1(ln.p) * (upper - lower) * (theta - lower) 4829 Pos <- (exp(ln.p) <= (upper- theta) / (upper - lower)) 4830 pstar <- (-expm1(ln.p) - (theta - lower) / (upper - lower)) / 4831 ((upper - theta) / (upper - lower)) 4832 } else { 4833 Neg <- (p >= (upper- theta) / (upper - lower)) 4834 temp1 <- (1 - p) * (upper - lower) * (theta - lower) 4835 Pos <- (p <= (upper- theta) / (upper - lower)) 4836 pstar <- ((upper- theta) / (upper - lower) - p) / 4837 ((upper - theta) / (upper - lower)) 4838 } 4839 } 4840 ans[ Neg] <- lower[ Neg] + sqrt(temp1[ Neg]) 4841 if (any(Pos)) { 4842 qstar <- cbind(1 - sqrt(1-pstar), 1 + sqrt(1-pstar)) 4843 qstar <- qstar[Pos,, drop = FALSE] 4844 qstar <- ifelse(qstar[, 1] >= 0 & qstar[, 1] <= 1, 4845 qstar[, 1], 4846 qstar[, 2]) 4847 ans[Pos] <- theta[Pos] + qstar * (upper - theta)[Pos] 4848 } 4849 4850 ans[theta < lower | theta > upper] <- NaN 4851 ans 4852} 4853 4854 4855 4856ptriangle <- function(q, theta, lower = 0, upper = 1, 4857 lower.tail = TRUE, log.p = FALSE) { 4858 4859 N <- max(length(q), length(theta), length(lower), length(upper)) 4860 if (length(q) != N) q <- rep_len(q, N) 4861 if (length(theta) != N) theta <- rep_len(theta, N) 4862 if (length(lower) != N) lower <- rep_len(lower, N) 4863 if (length(upper) != N) upper <- rep_len(upper, N) 4864 4865 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 4866 stop("bad input for argument 'lower.tail'") 4867 4868 if (!is.logical(log.p) || length(log.p) != 1) 4869 stop("bad input for argument 'log.p'") 4870 4871 ans <- q * 0 4872 qstar <- (q - lower)^2 / ((upper - lower) * (theta - lower)) 4873 Neg <- (lower <= q & q <= theta) 4874 4875 4876 ans[Neg] <- if (lower.tail) { 4877 if (log.p) { 4878 (log(qstar))[Neg] 4879 } else { 4880 qstar[Neg] 4881 } 4882 } else { 4883 if (log.p) { 4884 (log1p(-qstar))[Neg] 4885 } else { 4886 1 - qstar[Neg] 4887 } 4888 } 4889 4890 Pos <- (theta <= q & q <= upper) 4891 qstar <- (q - theta) / (upper-theta) 4892 4893 if (lower.tail) { 4894 if (log.p) { 4895 ans[Pos] <- log(((theta-lower)/(upper-lower))[Pos] + 4896 (qstar * (2-qstar) * 4897 (upper-theta) / (upper - lower))[Pos]) 4898 ans[q <= lower] <- -Inf 4899 ans[q >= upper] <- 0 4900 } else { 4901 ans[Pos] <- ((theta-lower)/(upper-lower))[Pos] + 4902 (qstar * (2-qstar) * 4903 (upper-theta) / (upper - lower))[Pos] 4904 ans[q <= lower] <- 0 4905 ans[q >= upper] <- 1 4906 } 4907 } else { 4908 if (log.p) { 4909 ans[Pos] <- log(((upper - theta)/(upper-lower))[Pos] + 4910 (qstar * (2-qstar) * 4911 (upper-theta) / (upper - lower))[Pos]) 4912 ans[q <= lower] <- 0 4913 ans[q >= upper] <- -Inf 4914 } else { 4915 ans[Pos] <- ((upper - theta)/(upper-lower))[Pos] + 4916 (qstar * (2-qstar) * 4917 (upper-theta) / (upper - lower))[Pos] 4918 ans[q <= lower] <- 1 4919 ans[q >= upper] <- 0 4920 } 4921 } 4922 4923 ans[theta < lower | theta > upper] <- NaN 4924 ans 4925} 4926 4927 4928 4929 4930 4931 4932 4933triangle.control <- function(stepsize = 0.33, maxit = 100, ...) { 4934 list(stepsize = stepsize, maxit = maxit) 4935} 4936 4937 4938 triangle <- 4939 function(lower = 0, upper = 1, 4940 link = extlogitlink(min = 0, max = 1), 4941 itheta = NULL) { 4942 4943 4944 4945 4946 4947 4948 if (!is.Numeric(lower)) 4949 stop("bad input for argument 'lower'") 4950 if (!is.Numeric(upper)) 4951 stop("bad input for argument 'upper'") 4952 if (!all(lower < upper)) 4953 stop("lower < upper values are required") 4954 4955 if (length(itheta) && !is.Numeric(itheta)) 4956 stop("bad input for 'itheta'") 4957 4958 4959 4960 4961 link <- as.list(substitute(link)) 4962 earg <- link2list(link) 4963 link <- attr(earg, "function.name") 4964 4965 4966 if (length(earg$min) && any(earg$min != lower)) 4967 stop("argument 'lower' does not match the 'link'") 4968 if (length(earg$max) && any(earg$max != upper)) 4969 stop("argument 'upper' does not match the 'link'") 4970 4971 4972 4973 new("vglmff", 4974 blurb = c("Triangle distribution\n\n", 4975 "Link: ", 4976 namesof("theta", link, earg = earg)), 4977 infos = eval(substitute(function(...) { 4978 list(M1 = 1, 4979 Q1 = 1, 4980 parameters.names = c("theta"), 4981 link = .link ) 4982 }, list( .link = link ))), 4983 4984 initialize = eval(substitute(expression({ 4985 4986 w.y.check(w = w, y = y, 4987 ncol.w.max = 1, 4988 ncol.y.max = 1) 4989 4990 4991 4992 4993 extra$lower <- rep_len( .lower , n) 4994 extra$upper <- rep_len( .upper , n) 4995 4996 if (any(y <= extra$lower | y >= extra$upper)) 4997 stop("some y values in [lower,upper] detected") 4998 4999 predictors.names <- 5000 namesof("theta", .link , earg = .earg , tag = FALSE) 5001 5002 5003 if (!length(etastart)) { 5004 Theta.init <- if (length( .itheta )) .itheta else { 5005 weighted.mean(y, w) 5006 } 5007 Theta.init <- rep_len(Theta.init, n) 5008 etastart <- theta2eta(Theta.init, .link , earg = .earg ) 5009 } 5010 }), list( .link = link, .earg = earg, .itheta=itheta, 5011 .upper = upper, .lower = lower ))), 5012 linkinv = eval(substitute(function(eta, extra = NULL) { 5013 Theta <- eta2theta(eta, .link , earg = .earg ) 5014 lower <- extra$lower 5015 upper <- extra$upper 5016 5017 mu1 <- (lower + upper + Theta) / 3 5018 5019 mu1 5020 }, list( .link = link, .earg = earg ))), 5021 last = eval(substitute(expression({ 5022 misc$link <- c(theta = .link ) 5023 5024 misc$earg <- list(theta = .earg ) 5025 5026 misc$expected <- TRUE 5027 }), list( .link = link, .earg = earg ))), 5028 loglikelihood = eval(substitute( 5029 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 5030 summation = TRUE) { 5031 Theta <- eta2theta(eta, .link , earg = .earg ) 5032 lower <- extra$lower 5033 upper <- extra$upper 5034 if (residuals) { 5035 stop("loglikelihood residuals not implemented yet") 5036 } else { 5037 ll.elts <- c(w) * dtriangle(x = y, theta = Theta, lower = lower, 5038 upper = upper, log = TRUE) 5039 if (summation) { 5040 sum(ll.elts) 5041 } else { 5042 ll.elts 5043 } 5044 } 5045 }, list( .link = link, .earg = earg ))), 5046 vfamily = c("triangle"), 5047 validparams = eval(substitute(function(eta, y, extra = NULL) { 5048 Theta <- eta2theta(eta, .link , earg = .earg ) 5049 okay1 <- all(is.finite(Theta)) && 5050 all(extra$lower < Theta & Theta < extra$upper) 5051 okay1 5052 }, list( .link = link, .earg = earg ))), 5053 5054 5055 5056 simslot = eval(substitute( 5057 function(object, nsim) { 5058 5059 pwts <- if (length(pwts <- object@prior.weights) > 0) 5060 pwts else weights(object, type = "prior") 5061 if (any(pwts != 1)) 5062 warning("ignoring prior weights") 5063 eta <- predict(object) 5064 extra <- object@extra 5065 Theta <- eta2theta(eta, .link , earg = .earg ) 5066 lower <- extra$lower 5067 upper <- extra$upper 5068 rtriangle(nsim * length(Theta), 5069 theta = Theta, lower = lower, upper = upper) 5070 }, list( .link = link, .earg = earg ))), 5071 5072 5073 5074 5075 deriv = eval(substitute(expression({ 5076 Theta <- eta2theta(eta, .link , earg = .earg ) 5077 dTheta.deta <- dtheta.deta(Theta, .link , earg = .earg ) 5078 5079 pos <- y > Theta 5080 neg <- y < Theta 5081 lower <- extra$lower 5082 upper <- extra$upper 5083 5084 dl.dTheta <- 0 * y 5085 dl.dTheta[neg] <- -1 / (Theta[neg]-lower[neg]) 5086 dl.dTheta[pos] <- 1 / (upper[pos]-Theta[pos]) 5087 5088 c(w) * dl.dTheta * dTheta.deta 5089 }), list( .link = link, .earg = earg ))), 5090 weight = eval(substitute(expression({ 5091 var.dl.dTheta <- 1 / ((Theta - lower) * (upper - Theta)) 5092 wz <- var.dl.dTheta * dTheta.deta^2 5093 c(w) * wz 5094 }), list( .link = link, .earg = earg )))) 5095} 5096 5097 5098 5099 5100 5101 5102 5103adjust0.loglaplace1 <- function(ymat, y, w, rep0) { 5104 rangey0 <- range(y[y > 0]) 5105 ymat[ymat <= 0] <- min(rangey0[1] / 2, rep0) 5106 ymat 5107} 5108 5109 5110loglaplace1.control <- function(maxit = 300, ...) { 5111 list(maxit = maxit) 5112} 5113 5114 5115 loglaplace1 <- function(tau = NULL, 5116 llocation = "loglink", 5117 ilocation = NULL, 5118 kappa = sqrt(tau/(1-tau)), 5119 Scale.arg = 1, 5120 ishrinkage = 0.95, 5121 parallel.locat = FALSE, digt = 4, 5122 idf.mu = 3, 5123 rep0 = 0.5, # 0.0001, 5124 minquantile = 0, maxquantile = Inf, 5125 imethod = 1, zero = NULL) { 5126 5127 if (length(minquantile) != 1) 5128 stop("bad input for argument 'minquantile'") 5129 if (length(maxquantile) != 1) 5130 stop("bad input for argument 'maxquantile'") 5131 5132 5133 if (!is.Numeric(rep0, positive = TRUE, length.arg = 1) || 5134 rep0 > 1) 5135 stop("bad input for argument 'rep0'") 5136 if (!is.Numeric(kappa, positive = TRUE)) 5137 stop("bad input for argument 'kappa'") 5138 5139 if (length(tau) && max(abs(kappa - sqrt(tau/(1-tau)))) > 1.0e-6) 5140 stop("arguments 'kappa' and 'tau' do not match") 5141 5142 5143 llocat <- as.list(substitute(llocation)) 5144 elocat <- link2list(llocat) 5145 llocat <- attr(elocat, "function.name") 5146 ilocat <- ilocation 5147 5148 5149 llocat.identity <- as.list(substitute("identitylink")) 5150 elocat.identity <- link2list(llocat.identity) 5151 llocat.identity <- attr(elocat.identity, "function.name") 5152 5153 5154 if (!is.Numeric(imethod, length.arg = 1, 5155 integer.valued = TRUE, positive = TRUE) || 5156 imethod > 4) 5157 stop("argument 'imethod' must be 1, 2 or ... 4") 5158 5159 5160 if (!is.Numeric(ishrinkage, length.arg = 1) || 5161 ishrinkage < 0 || 5162 ishrinkage > 1) 5163 stop("bad input for argument 'ishrinkage'") 5164 5165 if (!is.Numeric(Scale.arg, positive = TRUE)) 5166 stop("bad input for argument 'Scale.arg'") 5167 if (!is.logical(parallel.locat) || 5168 length(parallel.locat) != 1) 5169 stop("bad input for argument 'parallel.locat'") 5170 5171 fittedMean <- FALSE 5172 if (!is.logical(fittedMean) || length(fittedMean) != 1) 5173 stop("bad input for argument 'fittedMean'") 5174 5175 5176 mystring0 <- namesof("location", llocat, earg = elocat) 5177 mychars <- substring(mystring0, first = 1:nchar(mystring0), 5178 last = 1:nchar(mystring0)) 5179 mychars[nchar(mystring0)] <- ", inverse = TRUE)" 5180 mystring1 <- paste(mychars, collapse = "") 5181 5182 5183 5184 5185 new("vglmff", 5186 blurb = c("One-parameter ", 5187 if (llocat == "loglink") "log-Laplace" else 5188 c(llocat, "-Laplace"), 5189 " distribution\n\n", 5190 "Links: ", mystring0, "\n", "\n", 5191 "Quantiles: ", mystring1), 5192 constraints = eval(substitute(expression({ 5193 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 5194 bool = .parallel.locat , 5195 constraints = constraints, apply.int = FALSE) 5196 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 5197 predictors.names = predictors.names, 5198 M1 = 1) 5199 }), list( .parallel.locat = parallel.locat, 5200 .Scale.arg = Scale.arg, .zero = zero ))), 5201 5202 infos = eval(substitute(function(...) { 5203 list(M1 = 1, 5204 Q1 = 1, 5205 parameters.names = c("location"), 5206 llocation = .llocat ) 5207 }, list( .llocat = llocat, 5208 .zero = zero ))), 5209 5210 initialize = eval(substitute(expression({ 5211 extra$M <- M <- max(length( .Scale.arg ), length( .kappa )) 5212 extra$Scale <- rep_len( .Scale.arg , M) 5213 extra$kappa <- rep_len( .kappa , M) 5214 extra$tau <- extra$kappa^2 / (1 + extra$kappa^2) 5215 5216 5217 temp5 <- 5218 w.y.check(w = w, y = y, 5219 ncol.w.max = 1, 5220 ncol.y.max = 1, 5221 out.wy = TRUE, 5222 maximize = TRUE) 5223 w <- temp5$w 5224 y <- temp5$y 5225 5226 5227 5228 5229 extra$n <- n 5230 extra$y.names <- y.names <- 5231 paste("tau = ", round(extra$tau, digits = .digt), sep = "") 5232 extra$individual <- FALSE 5233 5234 5235 predictors.names <- 5236 namesof(paste("quantile(", y.names, ")", sep = ""), 5237 .llocat , earg = .elocat , tag = FALSE) 5238 5239 5240 if (FALSE) { 5241 if (min(y) < 0) 5242 stop("negative response values detected") 5243 if ((prop.0. <- weighted.mean(1*(y == 0), w)) >= min(extra$tau)) 5244 stop("sample proportion of 0s == ", round(prop.0., digits = 4), 5245 " > minimum 'tau' value. Choose larger values for 'tau'.") 5246 if ( .rep0 == 0.5 && 5247 (ave.tau <- (weighted.mean(1*(y <= 0), w) + 5248 weighted.mean(1*(y <= 1), w))/2) >= min(extra$tau)) 5249 warning("the minimum 'tau' value should be greater than ", 5250 round(ave.tau, digits = 4)) 5251 } 5252 5253 if (!length(etastart)) { 5254 if ( .imethod == 1) { 5255 locat.init <- quantile(rep(y, w), probs= extra$tau) + 1/16 5256 } else if ( .imethod == 2) { 5257 locat.init <- weighted.mean(y, w) 5258 } else if ( .imethod == 3) { 5259 locat.init <- median(y) 5260 } else if ( .imethod == 4) { 5261 Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 5262 y = y, w = w, 5263 df = .idf.mu ) 5264 locat.init <- c(predict(Fit5, x = x[, min(ncol(x), 2)])$y) 5265 } else { 5266 use.this <- weighted.mean(y, w) 5267 locat.init <- (1- .ishrinkage )*y + .ishrinkage * use.this 5268 } 5269 locat.init <- if (length( .ilocat )) 5270 rep_len( .ilocat , M) else 5271 rep_len(locat.init, M) 5272 locat.init <- matrix(locat.init, n, M, byrow = TRUE) 5273 if ( .llocat == "loglink") 5274 locat.init <- abs(locat.init) 5275 etastart <- 5276 cbind(theta2eta(locat.init, .llocat , earg = .elocat )) 5277 } 5278 }), list( .imethod = imethod, 5279 .idf.mu = idf.mu, .rep0 = rep0, 5280 .ishrinkage = ishrinkage, .digt = digt, 5281 .elocat = elocat, .Scale.arg = Scale.arg, 5282 .llocat = llocat, .kappa = kappa, 5283 .ilocat = ilocat ))), 5284 linkinv = eval(substitute(function(eta, extra = NULL) { 5285 locat.y = eta2theta(eta, .llocat , earg = .elocat ) 5286 if ( .fittedMean ) { 5287 stop("Yet to do: handle 'fittedMean = TRUE'") 5288 kappamat <- matrix(extra$kappa, extra$n, extra$M, byrow = TRUE) 5289 Scale <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 5290 locat.y + Scale * (1/kappamat - kappamat) 5291 } else { 5292 if (length(locat.y) > extra$n) 5293 dimnames(locat.y) <- list(dimnames(eta)[[1]], extra$y.names) 5294 locat.y 5295 } 5296 locat.y[locat.y < .minquantile] = .minquantile 5297 locat.y[locat.y > .maxquantile] = .maxquantile 5298 locat.y 5299 }, list( .elocat = elocat, .llocat = llocat, 5300 .minquantile = minquantile, .maxquantile = maxquantile, 5301 .fittedMean = fittedMean, .Scale.arg = Scale.arg, 5302 .kappa = kappa ))), 5303 last = eval(substitute(expression({ 5304 misc$link <- c(location = .llocat) 5305 5306 misc$earg <- list(location = .elocat ) 5307 5308 misc$expected <- TRUE 5309 5310 extra$kappa <- misc$kappa <- .kappa 5311 extra$tau <- misc$tau <- misc$kappa^2 / (1 + misc$kappa^2) 5312 extra$Scale.arg <- .Scale.arg 5313 5314 misc$true.mu <- .fittedMean # @fitted is not a true mu? 5315 misc$rep0 <- .rep0 5316 misc$minquantile <- .minquantile 5317 misc$maxquantile <- .maxquantile 5318 5319 extra$percentile <- numeric(length(misc$kappa)) 5320 locat.y <- as.matrix(locat.y) 5321 for (ii in seq_along(misc$kappa)) 5322 extra$percentile[ii] <- 100 * weighted.mean(y <= locat.y[, ii], w) 5323 }), list( .elocat = elocat, .llocat = llocat, 5324 .Scale.arg = Scale.arg, .fittedMean = fittedMean, 5325 .minquantile = minquantile, .maxquantile = maxquantile, 5326 .rep0 = rep0, .kappa = kappa ))), 5327 5328 loglikelihood = eval(substitute( 5329 function(mu, y, w, residuals = FALSE, eta, 5330 extra = NULL, 5331 summation = TRUE) { 5332 5333 kappamat <- matrix(extra$kappa, extra$n, extra$M, byrow = TRUE) 5334 Scale.w <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 5335 ymat <- matrix(y, extra$n, extra$M) 5336 5337 if ( .llocat == "loglink") 5338 ymat <- adjust0.loglaplace1(ymat = ymat, y = y, w = w, rep0 = .rep0) 5339 5340 w.mat <- theta2eta(ymat, .llocat , earg = .elocat ) # e.g., logofflink() 5341 5342 5343 5344 if (residuals) { 5345 stop("loglikelihood residuals not implemented yet") 5346 } else { 5347 ll.elts <- c(w) * dalap(x = c(w.mat), locat = c(eta), 5348 scale = c(Scale.w), kappa = c(kappamat), 5349 log = TRUE) 5350 if (summation) { 5351 sum(ll.elts) 5352 } else { 5353 ll.elts 5354 } 5355 } 5356 }, list( .elocat = elocat, .llocat = llocat, 5357 .rep0 = rep0, 5358 .Scale.arg = Scale.arg, .kappa = kappa ))), 5359 vfamily = c("loglaplace1"), 5360 validparams = eval(substitute(function(eta, y, extra = NULL) { 5361 locat.w <- eta 5362 locat.y <- eta2theta(locat.w, .llocat , earg = .elocat ) 5363 okay1 <- all(is.finite(locat.y)) 5364 okay1 5365 }, list( .elocat = elocat, .llocat = llocat, 5366 .rep0 = rep0, 5367 .Scale.arg = Scale.arg, .kappa = kappa ))), 5368 deriv = eval(substitute(expression({ 5369 ymat <- matrix(y, n, M) 5370 Scale.w <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 5371 locat.w <- eta 5372 locat.y <- eta2theta(locat.w, .llocat , earg = .elocat ) 5373 kappamat <- matrix(extra$kappa, n, M, byrow = TRUE) 5374 5375 ymat <- adjust0.loglaplace1(ymat = ymat, y = y, w = w, rep0= .rep0) 5376 w.mat <- theta2eta(ymat, .llocat , earg = .elocat ) # e.g., logitlink() 5377 zedd <- abs(w.mat-locat.w) / Scale.w 5378 dl.dlocat <- ifelse(w.mat >= locat.w, kappamat, 1/kappamat) * 5379 sqrt(2) * sign(w.mat-locat.w) / Scale.w 5380 5381 5382 dlocat.deta <- dtheta.deta(locat.w, 5383 .llocat.identity , 5384 earg = .elocat.identity ) 5385 c(w) * cbind(dl.dlocat * dlocat.deta) 5386 }), list( .Scale.arg = Scale.arg, .rep0 = rep0, 5387 .llocat = llocat, .elocat = elocat, 5388 .elocat.identity = elocat.identity, 5389 .llocat.identity = llocat.identity, 5390 5391 .kappa = kappa ))), 5392 weight = eval(substitute(expression({ 5393 ned2l.dlocat2 <- 2 / Scale.w^2 5394 wz <- cbind(ned2l.dlocat2 * dlocat.deta^2) 5395 c(w) * wz 5396 }), list( .Scale.arg = Scale.arg, 5397 .elocat = elocat, .llocat = llocat, 5398 .elocat.identity = elocat.identity, 5399 .llocat.identity = llocat.identity )))) 5400} 5401 5402 5403 5404 5405 5406loglaplace2.control <- function(save.weights = TRUE, ...) { 5407 list(save.weights = save.weights) 5408} 5409 5410 5411 loglaplace2 <- function(tau = NULL, 5412 llocation = "loglink", lscale = "loglink", 5413 ilocation = NULL, iscale = NULL, 5414 kappa = sqrt(tau/(1-tau)), 5415 ishrinkage = 0.95, 5416 parallel.locat = FALSE, digt = 4, 5417 eq.scale = TRUE, 5418 idf.mu = 3, 5419 rep0 = 0.5, nsimEIM = NULL, 5420 imethod = 1, zero = "(1 + M/2):M") { 5421 warning("it is best to use loglaplace1()") 5422 5423 if (length(nsimEIM) && 5424 (!is.Numeric(nsimEIM, length.arg = 1, 5425 integer.valued = TRUE) || 5426 nsimEIM <= 10)) 5427 stop("argument 'nsimEIM' should be an integer greater than 10") 5428 if (!is.Numeric(rep0, positive = TRUE, length.arg = 1) || 5429 rep0 > 1) 5430 stop("bad input for argument 'rep0'") 5431 if (!is.Numeric(kappa, positive = TRUE)) 5432 stop("bad input for argument 'kappa'") 5433 if (length(tau) && max(abs(kappa - sqrt(tau/(1-tau)))) > 1.0e-6) 5434 stop("arguments 'kappa' and 'tau' do not match") 5435 5436 5437 llocat <- as.list(substitute(llocation)) 5438 elocat <- link2list(llocat) 5439 llocat <- attr(elocat, "function.name") 5440 ilocat <- ilocation 5441 5442 lscale <- as.list(substitute(lscale)) 5443 escale <- link2list(lscale) 5444 lscale <- attr(escale, "function.name") 5445 5446 5447 5448 5449 if (!is.Numeric(imethod, length.arg = 1, 5450 integer.valued = TRUE, positive = TRUE) || 5451 imethod > 4) 5452 stop("argument 'imethod' must be 1, 2 or ... 4") 5453 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 5454 stop("bad input for argument 'iscale'") 5455 5456 5457 if (!is.Numeric(ishrinkage, length.arg = 1) || 5458 ishrinkage < 0 || 5459 ishrinkage > 1) 5460 stop("bad input for argument 'ishrinkage'") 5461 if (!is.logical(eq.scale) || length(eq.scale) != 1) 5462 stop("bad input for argument 'eq.scale'") 5463 if (!is.logical(parallel.locat) || 5464 length(parallel.locat) != 1) 5465 stop("bad input for argument 'parallel.locat'") 5466 fittedMean <- FALSE 5467 if (!is.logical(fittedMean) || length(fittedMean) != 1) 5468 stop("bad input for argument 'fittedMean'") 5469 5470 if (llocat != "loglink") 5471 stop("argument 'llocat' must be \"loglink\"") 5472 5473 5474 new("vglmff", 5475 blurb = c("Two-parameter log-Laplace distribution\n\n", 5476 "Links: ", 5477 namesof("location", llocat, earg = elocat), ", ", 5478 namesof("scale", lscale, earg = escale), 5479 "\n", "\n", 5480 "Mean: zz location + scale * ", 5481 "(1/kappa - kappa) / sqrt(2)", "\n", 5482 "Quantiles: location", "\n", 5483 "Variance: zz scale^2 * (1 + kappa^4) / (2 * kappa^2)"), 5484 constraints = eval(substitute(expression({ 5485 .ZERO <- .zero 5486 if (is.character( .ZERO )) 5487 .ZERO <- eval(parse(text = .ZERO )) 5488 .PARALLEL <- .parallel.locat 5489 parelHmat <- if (is.logical( .PARALLEL ) && .PARALLEL ) 5490 matrix(1, M/2, 1) else diag(M/2) 5491 scaleHmat <- if (is.logical( .eq.scale ) && .eq.scale ) 5492 matrix(1, M/2, 1) else diag(M/2) 5493 mycmatrix <- cbind(rbind( parelHmat, 0*parelHmat), 5494 rbind(0*scaleHmat, scaleHmat)) 5495 constraints <- cm.VGAM(mycmatrix, x = x, 5496 bool = .PARALLEL , 5497 constraints = constraints, 5498 apply.int = FALSE) 5499 constraints <- cm.zero.VGAM(constraints, x = x, .ZERO , M = M, 5500 predictors.names = predictors.names, 5501 M1 = 2) 5502 5503 if ( .PARALLEL && names(constraints)[1] == "(Intercept)") { 5504 parelHmat <- diag(M/2) 5505 mycmatrix <- cbind(rbind( parelHmat, 0*parelHmat), 5506 rbind(0*scaleHmat, scaleHmat)) 5507 constraints[["(Intercept)"]] <- mycmatrix 5508 } 5509 if (is.logical( .eq.scale) && .eq.scale && 5510 names(constraints)[1] == "(Intercept)") { 5511 temp3 <- constraints[["(Intercept)"]] 5512 temp3 <- cbind(temp3[,1:(M/2)], rbind(0*scaleHmat, scaleHmat)) 5513 constraints[["(Intercept)"]] = temp3 5514 } 5515 }), list( .eq.scale = eq.scale, .parallel.locat = parallel.locat, 5516 .zero = zero ))), 5517 initialize = eval(substitute(expression({ 5518 extra$kappa <- .kappa 5519 extra$tau <- extra$kappa^2 / (1 + extra$kappa^2) 5520 5521 5522 5523 temp5 <- 5524 w.y.check(w = w, y = y, 5525 ncol.w.max = 1, 5526 ncol.y.max = 1, 5527 out.wy = TRUE, 5528 maximize = TRUE) 5529 w <- temp5$w 5530 y <- temp5$y 5531 5532 5533 5534 5535 extra$M <- M <- 2 * length(extra$kappa) 5536 extra$n <- n 5537 extra$y.names <- y.names <- 5538 paste("tau = ", round(extra$tau, digits = .digt), sep = "") 5539 extra$individual = FALSE 5540 5541 predictors.names <- 5542 c(namesof(paste("quantile(", y.names, ")", sep = ""), 5543 .llocat , earg = .elocat, tag = FALSE), 5544 namesof(if (M == 2) "scale" else 5545 paste("scale", 1:(M/2), sep = ""), 5546 .lscale , earg = .escale, tag = FALSE)) 5547 if (weighted.mean(1 * (y < 0.001), w) >= min(extra$tau)) 5548 stop("sample proportion of 0s > minimum 'tau' value. ", 5549 "Choose larger values for 'tau'.") 5550 5551 if (!length(etastart)) { 5552 if ( .imethod == 1) { 5553 locat.init.y <- weighted.mean(y, w) 5554 scale.init <- sqrt(var(y) / 2) 5555 } else if ( .imethod == 2) { 5556 locat.init.y <- median(y) 5557 scale.init <- sqrt(sum(c(w)*abs(y-median(y))) / (sum(w) *2)) 5558 } else if ( .imethod == 3) { 5559 Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)], 5560 y = y, w = w, 5561 df = .idf.mu ) 5562 locat.init.y <- c(predict(Fit5, x = x[, min(ncol(x), 2)])$y) 5563 scale.init <- sqrt(sum(c(w)*abs(y-median(y))) / (sum(w) *2)) 5564 } else { 5565 use.this <- weighted.mean(y, w) 5566 locat.init.y <- (1- .ishrinkage )*y + .ishrinkage * use.this 5567 scale.init <- sqrt(sum(c(w)*abs(y-median(y ))) / (sum(w) *2)) 5568 } 5569 locat.init.y <- if (length( .ilocat )) 5570 rep_len( .ilocat , n) else 5571 rep_len(locat.init.y, n) 5572 locat.init.y <- matrix(locat.init.y, n, M/2) 5573 scale.init <- if (length( .iscale )) 5574 rep_len( .iscale , n) else 5575 rep_len(scale.init, n) 5576 scale.init <- matrix(scale.init, n, M/2) 5577 etastart <- 5578 cbind(theta2eta(locat.init.y, .llocat , earg = .elocat ), 5579 theta2eta(scale.init, .lscale , earg = .escale )) 5580 } 5581 }), list( .imethod = imethod, 5582 .idf.mu = idf.mu, .kappa = kappa, 5583 .ishrinkage = ishrinkage, .digt = digt, 5584 .llocat = llocat, .lscale = lscale, 5585 .elocat = elocat, .escale = escale, 5586 .ilocat = ilocat, .iscale = iscale ))), 5587 linkinv = eval(substitute(function(eta, extra = NULL) { 5588 locat.y <- eta2theta(eta[, 1:(extra$M/2), drop = FALSE], 5589 .llocat , earg = .elocat ) 5590 if ( .fittedMean ) { 5591 kappamat <- matrix(extra$kappa, extra$n, extra$M/2, 5592 byrow = TRUE) 5593 Scale.y <- eta2theta(eta[,(1+extra$M/2):extra$M], 5594 .lscale , earg = .escale ) 5595 locat.y + Scale.y * (1/kappamat - kappamat) 5596 } else { 5597 dimnames(locat.y) = list(dimnames(eta)[[1]], extra$y.names) 5598 locat.y 5599 } 5600 }, list( .llocat = llocat, .lscale = lscale, 5601 .elocat = elocat, .escale = escale, 5602 .fittedMean = fittedMean, 5603 .kappa = kappa ))), 5604 last = eval(substitute(expression({ 5605 misc$link <- c(location = .llocat , scale = .lscale ) 5606 5607 misc$earg <- list(location = .elocat , scale = .escale ) 5608 5609 misc$expected <- TRUE 5610 extra$kappa <- misc$kappa <- .kappa 5611 extra$tau <- misc$tau <- misc$kappa^2 / (1 + misc$kappa^2) 5612 misc$true.mu <- .fittedMean # @fitted is not a true mu? 5613 misc$nsimEIM <- .nsimEIM 5614 misc$rep0 <- .rep0 5615 extra$percentile <- numeric(length(misc$kappa)) 5616 locat <- as.matrix(locat.y) 5617 for (ii in seq_along(misc$kappa)) 5618 extra$percentile[ii] <- 100 * 5619 weighted.mean(y <= locat.y[, ii], w) 5620 }), list( .elocat = elocat, .llocat = llocat, 5621 .escale = escale, .lscale = lscale, 5622 .fittedMean = fittedMean, 5623 .nsimEIM = nsimEIM, .rep0 = rep0, 5624 .kappa = kappa ))), 5625 5626 loglikelihood = eval(substitute( 5627 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 5628 summation = TRUE) { 5629 5630 kappamat <- matrix(extra$kappa, extra$n, extra$M/2, byrow = TRUE) 5631 Scale.w <- eta2theta(eta[, (1+extra$M/2):extra$M], 5632 .lscale , earg = .escale ) 5633 ymat <- matrix(y, extra$n, extra$M/2) 5634 ymat[ymat <= 0] <- min(min(y[y > 0]), .rep0 ) # Adjust for 0s 5635 ell.mat <- matrix(c(dloglaplace(x = c(ymat), 5636 locat.ald = c(eta[, 1:(extra$M/2)]), 5637 scale.ald = c(Scale.w), 5638 kappa = c(kappamat), log = TRUE)), 5639 extra$n, extra$M/2) 5640 if (residuals) { 5641 stop("loglikelihood residuals not implemented yet") 5642 } else { 5643 ll.elts <- c(w) * ell.mat 5644 if (summation) { 5645 sum(ll.elts) 5646 } else { 5647 ll.elts 5648 } 5649 } 5650 }, list( .elocat = elocat, .llocat = llocat, 5651 .escale = escale, .lscale = lscale, 5652 .rep0 = rep0, .kappa = kappa ))), 5653 5654 5655 vfamily = c("loglaplace2"), 5656 validparams = eval(substitute(function(eta, y, extra = NULL) { 5657 Scale.w <- eta2theta(eta[, (1+extra$M/2):extra$M], 5658 .lscale , earg = .escale ) 5659 locat.w <- eta[, 1:(extra$M/2), drop = FALSE] 5660 locat.y <- eta2theta(locat.w, .llocat , earg = .elocat ) 5661 okay1 <- all(is.finite(locat.y)) && 5662 all(is.finite(Scale.w)) && all(0 < Scale.w) 5663 okay1 5664 }, list( .elocat = elocat, .llocat = llocat, 5665 .escale = escale, .lscale = lscale, 5666 .rep0 = rep0, .kappa = kappa ))), 5667 deriv = eval(substitute(expression({ 5668 ymat <- matrix(y, n, M/2) 5669 Scale.w <- eta2theta(eta[, (1+extra$M/2):extra$M], 5670 .lscale , earg = .escale ) 5671 locat.w <- eta[, 1:(extra$M/2), drop = FALSE] 5672 locat.y <- eta2theta(locat.w, .llocat , earg = .elocat ) 5673 kappamat <- matrix(extra$kappa, n, M/2, byrow = TRUE) 5674 w.mat <- ymat 5675 w.mat[w.mat <= 0] <- min(min(w.mat[w.mat > 0]), .rep0) 5676 w.mat <- theta2eta(w.mat, .llocat , earg = .elocat ) 5677 zedd <- abs(w.mat-locat.w) / Scale.w 5678 dl.dlocat <- sqrt(2) * 5679 ifelse(w.mat >= locat.w, kappamat, 1/kappamat) * 5680 sign(w.mat-locat.w) / Scale.w 5681 dl.dscale <- sqrt(2) * 5682 ifelse(w.mat >= locat.w, kappamat, 1/kappamat) * 5683 zedd / Scale.w - 1 / Scale.w 5684 dlocat.deta <- dtheta.deta(locat.w, .llocat , earg = .elocat ) 5685 dscale.deta <- dtheta.deta(Scale.w, .lscale , earg = .escale ) 5686 c(w) * cbind(dl.dlocat * dlocat.deta, 5687 dl.dscale * dscale.deta) 5688 }), list( .escale = escale, .lscale = lscale, 5689 .elocat = elocat, .llocat = llocat, 5690 .rep0 = rep0, .kappa = kappa ))), 5691 weight = eval(substitute(expression({ 5692 run.varcov <- 0 5693 ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE) 5694 dthetas.detas <- cbind(dlocat.deta, dscale.deta) 5695 if (length( .nsimEIM )) { 5696 for (ii in 1:( .nsimEIM )) { 5697 wsim <- matrix(rloglap(n*M/2, loc = c(locat.w), 5698 sca = c(Scale.w), 5699 kappa = c(kappamat)), n, M/2) 5700 zedd <- abs(wsim-locat.w) / Scale.w 5701 dl.dlocat <- sqrt(2) * 5702 ifelse(wsim >= locat.w, kappamat, 1/kappamat) * 5703 sign(wsim-locat.w) / Scale.w 5704 dl.dscale <- sqrt(2) * 5705 ifelse(wsim >= locat.w, kappamat, 1/kappamat) * 5706 zedd / Scale.w - 1 / Scale.w 5707 5708 rm(wsim) 5709 temp3 <- cbind(dl.dlocat, dl.dscale) # n x M matrix 5710 run.varcov <- ((ii-1) * run.varcov + 5711 temp3[,ind1$row.index]*temp3[,ind1$col.index]) / ii 5712 } 5713 wz <- if (intercept.only) 5714 matrix(colMeans(run.varcov), 5715 n, ncol(run.varcov), byrow = TRUE) else run.varcov 5716 5717 wz <- wz * dthetas.detas[,ind1$row] * dthetas.detas[,ind1$col] 5718 wz <- c(w) * matrix(wz, n, dimm(M)) 5719 wz 5720 } else { 5721 d2l.dlocat2 <- 2 / (Scale.w * locat.w)^2 5722 d2l.dscale2 <- 1 / Scale.w^2 5723 wz <- cbind(d2l.dlocat2 * dlocat.deta^2, 5724 d2l.dscale2 * dscale.deta^2) 5725 c(w) * wz 5726 } 5727 }), list( .elocat = elocat, .escale = escale, 5728 .llocat = llocat, .lscale = lscale, 5729 .nsimEIM = nsimEIM) ))) 5730} 5731 5732 5733 5734 5735 5736 5737logitlaplace1.control <- function(maxit = 300, ...) { 5738 list(maxit = maxit) 5739} 5740 5741 5742adjust01.logitlaplace1 <- function(ymat, y, w, rep01) { 5743 rangey01 <- range(y[(y > 0) & (y < 1)]) 5744 ymat[ymat <= 0] <- min(rangey01[1] / 2, rep01 / w[y <= 0]) 5745 ymat[ymat >= 1] <- max((1 + rangey01[2]) / 2, 1 - rep01 / w[y >= 1]) 5746 ymat 5747} 5748 5749 5750 5751 5752 5753 logitlaplace1 <- 5754 function(tau = NULL, 5755 llocation = "logitlink", 5756 ilocation = NULL, 5757 kappa = sqrt(tau/(1-tau)), 5758 Scale.arg = 1, 5759 ishrinkage = 0.95, parallel.locat = FALSE, digt = 4, 5760 idf.mu = 3, 5761 rep01 = 0.5, 5762 imethod = 1, zero = NULL) { 5763 5764 if (!is.Numeric(rep01, positive = TRUE, length.arg = 1) || 5765 rep01 > 0.5) 5766 stop("bad input for argument 'rep01'") 5767 if (!is.Numeric(kappa, positive = TRUE)) 5768 stop("bad input for argument 'kappa'") 5769 5770 if (length(tau) && max(abs(kappa - sqrt(tau/(1-tau)))) > 1.0e-6) 5771 stop("arguments 'kappa' and 'tau' do not match") 5772 5773 5774 llocat <- as.list(substitute(llocation)) 5775 elocat <- link2list(llocat) 5776 llocat <- attr(elocat, "function.name") 5777 ilocat <- ilocation 5778 5779 5780 llocat.identity <- as.list(substitute("identitylink")) 5781 elocat.identity <- link2list(llocat.identity) 5782 llocat.identity <- attr(elocat.identity, "function.name") 5783 5784 5785 5786 5787 if (!is.Numeric(imethod, length.arg = 1, 5788 integer.valued = TRUE, positive = TRUE) || 5789 imethod > 4) 5790 stop("argument 'imethod' must be 1, 2 or ... 4") 5791 5792 if (!is.Numeric(ishrinkage, length.arg = 1) || 5793 ishrinkage < 0 || 5794 ishrinkage > 1) 5795 stop("bad input for argument 'ishrinkage'") 5796 5797 if (!is.Numeric(Scale.arg, positive = TRUE)) 5798 stop("bad input for argument 'Scale.arg'") 5799 if (!is.logical(parallel.locat) || 5800 length(parallel.locat) != 1) 5801 stop("bad input for argument 'parallel.locat'") 5802 fittedMean <- FALSE 5803 if (!is.logical(fittedMean) || 5804 length(fittedMean) != 1) 5805 stop("bad input for argument 'fittedMean'") 5806 5807 5808 mystring0 <- namesof("location", llocat, earg = elocat) 5809 mychars <- substring(mystring0, first = 1:nchar(mystring0), 5810 last = 1:nchar(mystring0)) 5811 mychars[nchar(mystring0)] = ", inverse = TRUE)" 5812 mystring1 <- paste(mychars, collapse = "") 5813 5814 5815 5816 5817 new("vglmff", 5818 blurb = c("One-parameter ", llocat, "-Laplace distribution\n\n", 5819 "Links: ", mystring0, "\n", "\n", 5820 "Quantiles: ", mystring1), 5821 constraints = eval(substitute(expression({ 5822 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 5823 bool = .parallel.locat , 5824 constraints = constraints, apply.int = FALSE) 5825 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 5826 predictors.names = predictors.names, 5827 M1 = 1) 5828 }), list( .parallel.locat = parallel.locat, 5829 .Scale.arg = Scale.arg, .zero = zero ))), 5830 5831 infos = eval(substitute(function(...) { 5832 list(M1 = 1, 5833 Q1 = 1, 5834 multipleResponses = FALSE, 5835 parameters.names = c("location"), 5836 llocation = .llocat , 5837 zero = .zero ) 5838 }, list( .zero = zero, 5839 .llocat = llocat ))), 5840 5841 initialize = eval(substitute(expression({ 5842 extra$M <- M <- max(length( .Scale.arg ), length( .kappa )) 5843 extra$Scale <- rep_len( .Scale.arg , M) 5844 extra$kappa <- rep_len( .kappa , M) 5845 extra$tau <- extra$kappa^2 / (1 + extra$kappa^2) 5846 5847 5848 5849 temp5 <- 5850 w.y.check(w = w, y = y, 5851 ncol.w.max = 1, 5852 ncol.y.max = 1, 5853 out.wy = TRUE, 5854 maximize = TRUE) 5855 w <- temp5$w 5856 y <- temp5$y 5857 5858 5859 5860 5861 5862 extra$n <- n 5863 extra$y.names <- y.names <- 5864 paste("tau = ", round(extra$tau, digits = .digt), sep = "") 5865 extra$individual <- FALSE 5866 5867 predictors.names <- 5868 namesof(paste("quantile(", y.names, ")", sep = ""), 5869 .llocat , earg = .elocat, tag = FALSE) 5870 5871 if (all(y == 0 | y == 1)) 5872 stop("response cannot be all 0s or 1s") 5873 if (min(y) < 0) 5874 stop("negative response values detected") 5875 if (max(y) > 1) 5876 stop("response values greater than 1 detected") 5877 if ((prop.0. <- weighted.mean(1*(y == 0), w)) >= min(extra$tau)) 5878 stop("sample proportion of 0s == ", round(prop.0., digits = 4), 5879 " > minimum 'tau' value. Choose larger values for 'tau'.") 5880 if ((prop.1. <- weighted.mean(1*(y == 1), w)) >= max(extra$tau)) 5881 stop("sample proportion of 1s == ", round(prop.1., digits = 4), 5882 " < maximum 'tau' value. Choose smaller values for 'tau'.") 5883 if (!length(etastart)) { 5884 if ( .imethod == 1) { 5885 locat.init <- quantile(rep(y, w), probs= extra$tau) 5886 } else if ( .imethod == 2) { 5887 locat.init <- weighted.mean(y, w) 5888 locat.init <- median(rep(y, w)) 5889 } else if ( .imethod == 3) { 5890 use.this <- weighted.mean(y, w) 5891 locat.init <- (1- .ishrinkage )*y + use.this * .ishrinkage 5892 } else { 5893 stop("this option not implemented") 5894 } 5895 5896 5897 locat.init <- if (length( .ilocat )) 5898 rep_len( .ilocat , M) else 5899 rep_len(locat.init, M) 5900 locat.init <- matrix(locat.init, n, M, byrow = TRUE) 5901 locat.init <- abs(locat.init) 5902 etastart <- 5903 cbind(theta2eta(locat.init, .llocat , earg = .elocat )) 5904 } 5905 }), list( .imethod = imethod, 5906 .idf.mu = idf.mu, 5907 .ishrinkage = ishrinkage, .digt = digt, 5908 .elocat = elocat, .Scale.arg = Scale.arg, 5909 .llocat = llocat, .kappa = kappa, 5910 .ilocat = ilocat ))), 5911 linkinv = eval(substitute(function(eta, extra = NULL) { 5912 locat.y <- eta2theta(eta, .llocat , earg = .elocat ) 5913 if ( .fittedMean ) { 5914 stop("Yet to do: handle 'fittedMean = TRUE'") 5915 kappamat <- matrix(extra$kappa, extra$n, extra$M, byrow = TRUE) 5916 Scale <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 5917 locat.y + Scale * (1/kappamat - kappamat) 5918 } else { 5919 if (length(locat.y) > extra$n) 5920 dimnames(locat.y) <- list(dimnames(eta)[[1]], extra$y.names) 5921 locat.y 5922 } 5923 }, list( .elocat = elocat, .llocat = llocat, 5924 .fittedMean = fittedMean, .Scale.arg = Scale.arg, 5925 .kappa = kappa ))), 5926 last = eval(substitute(expression({ 5927 misc$link <- c(location = .llocat ) 5928 misc$earg <- list(location = .elocat ) 5929 5930 misc$expected <- TRUE 5931 5932 extra$kappa <- misc$kappa <- .kappa 5933 extra$tau <- misc$tau <- misc$kappa^2 / (1 + misc$kappa^2) 5934 extra$Scale.arg <- .Scale.arg 5935 5936 misc$true.mu <- .fittedMean # @fitted is not a true mu? 5937 misc$rep01 <- .rep01 5938 5939 extra$percentile <- numeric(length(misc$kappa)) 5940 locat.y <- eta2theta(eta, .llocat , earg = .elocat ) 5941 locat.y <- as.matrix(locat.y) 5942 for (ii in seq_along(misc$kappa)) 5943 extra$percentile[ii] <- 100 * 5944 weighted.mean(y <= locat.y[, ii], w) 5945 5946 }), list( .elocat = elocat, .llocat = llocat, 5947 .Scale.arg = Scale.arg, .fittedMean = fittedMean, 5948 .rep01 = rep01, 5949 .kappa = kappa ))), 5950 5951 5952 loglikelihood = eval(substitute( 5953 function(mu, y, w, residuals = FALSE, eta, 5954 extra = NULL, 5955 summation = TRUE) { 5956 5957 kappamat <- matrix(extra$kappa, extra$n, extra$M, byrow = TRUE) 5958 Scale.w <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 5959 ymat <- matrix(y, extra$n, extra$M) 5960 ymat <- adjust01.logitlaplace1(ymat = ymat, y = y, w = w, 5961 rep01 = .rep01) 5962 w.mat <- theta2eta(ymat, .llocat , earg = .elocat ) # e.g., logitlink() 5963 if (residuals) { 5964 stop("loglikelihood residuals not implemented yet") 5965 } else { 5966 ll.elts <- 5967 c(w) * dalap(x = c(w.mat), location = c(eta), 5968 scale = c(Scale.w), kappa = c(kappamat), 5969 log = TRUE) 5970 if (summation) { 5971 sum(ll.elts) 5972 } else { 5973 ll.elts 5974 } 5975 } 5976 }, list( .elocat = elocat, .llocat = llocat, 5977 .rep01 = rep01, 5978 .Scale.arg = Scale.arg, .kappa = kappa ))), 5979 5980 5981 vfamily = c("logitlaplace1"), 5982 validparams = eval(substitute(function(eta, y, extra = NULL) { 5983 locat.w <- eta 5984 okay1 <- all(is.finite(locat.w)) 5985 okay1 5986 }, list( .Scale.arg = Scale.arg, .rep01 = rep01, 5987 .elocat = elocat, 5988 .llocat = llocat, 5989 5990 .elocat.identity = elocat.identity, 5991 .llocat.identity = llocat.identity, 5992 5993 .kappa = kappa ))), 5994 deriv = eval(substitute(expression({ 5995 ymat <- matrix(y, n, M) 5996 Scale.w <- matrix(extra$Scale, extra$n, extra$M, byrow = TRUE) 5997 locat.w <- eta 5998 kappamat <- matrix(extra$kappa, n, M, byrow = TRUE) 5999 ymat <- adjust01.logitlaplace1(ymat = ymat, y = y, w = w, 6000 rep01 = .rep01 ) 6001 w.mat <- theta2eta(ymat, .llocat , earg = .elocat ) # e.g., logitlink() 6002 zedd <- abs(w.mat - locat.w) / Scale.w 6003 dl.dlocat <- ifelse(w.mat >= locat.w, kappamat, 1/kappamat) * 6004 sqrt(2) * sign(w.mat-locat.w) / Scale.w 6005 6006 6007 dlocat.deta <- dtheta.deta(locat.w, 6008 "identitylink", 6009 earg = .elocat.identity ) 6010 6011 6012 c(w) * cbind(dl.dlocat * dlocat.deta) 6013 }), list( .Scale.arg = Scale.arg, .rep01 = rep01, 6014 .elocat = elocat, 6015 .llocat = llocat, 6016 6017 .elocat.identity = elocat.identity, 6018 .llocat.identity = llocat.identity, 6019 6020 .kappa = kappa ))), 6021 weight = eval(substitute(expression({ 6022 d2l.dlocat2 <- 2 / Scale.w^2 6023 wz <- cbind(d2l.dlocat2 * dlocat.deta^2) 6024 c(w) * wz 6025 }), list( .Scale.arg = Scale.arg, 6026 .elocat = elocat, .llocat = llocat )))) 6027} 6028 6029 6030 6031 6032 6033 6034 6035 6036 6037 6038 6039 6040 6041 6042 6043 6044 6045 6046 6047 6048 6049