1# This program is free software; you can redistribute it and/or modify 2# it under the terms of the GNU General Public License as published by 3# the Free Software Foundation; either version 2 of the License, or 4# (at your option) any later version. 5# 6# This program is distributed in the hope that it will be useful, 7# but WITHOUT ANY WARRANTY; without even the implied warranty of 8# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 9# GNU General Public License for more details. 10# 11# A copy of the GNU General Public License is available at 12# http://www.r-project.org/Licenses/ 13 14 15.invTest <- function(object, which, level = 0.95, fact = 3, type=c("LR", "LM", "J"), corr=NULL) 16 { 17 type <- match.arg(type) 18 argCall <- object$allArg 19 if (length(which) > 1) 20 stop("tests are inverted only for one parameter") 21 if (is.character(which)) 22 { 23 which <- which(which == names(object$coefficients)) 24 if (length(which) == 0) 25 stop("Wrong parameter names") 26 } 27 wtest <- which(type==c("LR", "LM", "J")) 28 if (object$df > 0) 29 { 30 test0 <- specTest(object)$test[wtest,] 31 test0 <- test0[1] 32 } else { 33 test0 <- 0 34 } 35 if (length(object$coefficients) == 1) 36 { 37 argCall$optfct <- NULL 38 argCall$constraint <- NULL 39 argCall$eqConst <- NULL 40 argCall$eqConstFullVcov <- NULL 41 fctCall <- "evalGel" 42 } else { 43 fctCall <- "gel" 44 argCall$eqConst <- which 45 } 46 if (length(object$coefficients) == 2) 47 { 48 sdcoef <- sqrt(diag(vcov(object))[-which]) 49 coef <- object$coefficients[-which] 50 upper <- coef + fact*sdcoef 51 lower <- coef - fact*sdcoef 52 argCall$method <- "Brent" 53 argCall$lower <- lower 54 argCall$upper <- upper 55 } 56 sdcoef <- sqrt(diag(vcov(object))[which]) 57 coef <- object$coefficients[which] 58 int1 <- c(coef, coef + fact*sdcoef) 59 int2 <- c(coef - fact*sdcoef, coef) 60 fct <- function(coef, which, wtest, level, test0, corr=NULL) 61 { 62 argCall$tet0 <- object$coefficients 63 argCall$tet0[which] <- coef 64 obj <- do.call(get(fctCall), argCall) 65 test <- as.numeric(specTest(obj)$test[wtest,1]) - test0 66 if (is.null(corr)) 67 level - pchisq(test, 1) 68 else 69 level - pchisq(test/corr, 1) 70 } 71 res1 <- uniroot(fct, int1, which = which, wtest=wtest, level=level, 72 test0=test0, corr=corr) 73 res2 <- uniroot(fct, int2, which = which, wtest=wtest, level=level, 74 test0=test0, corr=corr) 75 sort(c(res1$root, res2$root)) 76 } 77 78 79confint.gel <- function(object, parm, level = 0.95, lambda = FALSE, 80 type = c("Wald", "invLR", "invLM", "invJ"), 81 fact = 3, corr = NULL, ...) 82 { 83 type <- match.arg(type) 84 z <- object 85 n <- nrow(z$gt) 86 if (type == "Wald") 87 { 88 ntest <- "Direct Wald type confidence interval" 89 se_par <- sqrt(diag(z$vcov_par)) 90 par <- z$coefficients 91 tval <- par/se_par 92 se_parl <- sqrt(diag(z$vcov_lambda)) 93 lamb <- z$lambda 94 zs <- qnorm((1 - level)/2, lower.tail=FALSE) 95 ch <- zs*se_par 96 if(!lambda) 97 { 98 ans <- cbind(par-ch, par+ch) 99 dimnames(ans) <- list(names(par), c((1 - level)/2, 0.5+level/2)) 100 } 101 if(lambda) 102 { 103 if (length(z$coefficients) == length(z$lambda)) 104 { 105 cat("\nNo confidence intervals for lambda when the model is just identified.\n") 106 return(NULL) 107 } else { 108 chl <- zs*se_parl 109 ans <- cbind(lamb - chl, lamb + chl) 110 dimnames(ans) <- list(names(lamb), c((1 - level)/2, 0.5 + level/2)) 111 } 112 } 113 if(!missing(parm)) 114 ans <- ans[parm,] 115 } else { 116 if(missing(parm)) 117 parm <- names(object$coefficients) 118 type <- strsplit(type, "v")[[1]][2] 119 ntest <- paste("Confidence interval based on the inversion of the ", type, " test", sep="") 120 ans <- lapply(parm, function(w) .invTest(object, w, level = level, fact = fact, type=type, corr=corr)) 121 ans <- do.call(rbind, ans) 122 if (!is.character(parm)) 123 parm <- names(object$coefficients)[parm] 124 dimnames(ans) <- list(parm, c((1 - level)/2, 0.5+level/2)) 125 } 126 ans <- list(test=ans,ntest=ntest) 127 class(ans) <- "confint" 128 ans 129 } 130 131print.confint <- function(x, digits = 5, ...) 132 { 133 cat("\n", x$ntest, sep="") 134 cat("\n#######################################\n") 135 print.default(format(x$test, digits = digits), 136 print.gap = 2, quote = FALSE) 137 invisible(x) 138 } 139 140coef.gel <- function(object, lambda = FALSE, ...) 141 { 142 if(!lambda) 143 object$coefficients 144 else 145 object$lambda 146 } 147 148vcov.gel <- function(object, lambda = FALSE, ...) 149 { 150 if(!lambda) 151 object$vcov_par 152 else 153 object$vcov_lambda 154 } 155 156print.gel <- function(x, digits = 5, ...) 157 { 158 if (is.null(x$CGEL)) 159 cat("Type de GEL: ", x$typeDesc, "\n") 160 else 161 cat("CGEL of type: ", x$typeDesc, " (alpha = ", x$CGEL, ")\n") 162 if (!is.null(attr(x$dat,"smooth"))) 163 { 164 cat("Kernel: ", attr(x$dat,"smooth")$kernel," (bw=", 165 attr(x$dat,"smooth")$bw,")\n\n") 166 } 167 else 168 cat("\n") 169 170 cat("Coefficients:\n") 171 print.default(format(coef(x), digits = digits), 172 print.gap = 2, quote = FALSE) 173 cat("\n") 174 cat("Lambdas:\n") 175 print.default(format(coef(x, lambda = TRUE), digits = digits), 176 print.gap = 2, quote = FALSE) 177 cat("\n") 178 cat("Convergence code for the coefficients: ", x$conv_par,"\n") 179 cat("Convergence code for Lambda: ", x$conv_lambda$convergence,"\n") 180 cat(x$specMod) 181 invisible(x) 182 } 183 184print.summary.gel <- function(x, digits = 5, ...) 185 { 186 cat("\nCall:\n") 187 cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n", sep = "") 188 if (is.null(x$CGEL)) 189 cat("Type of GEL: ", x$typeDesc, "\n") 190 else 191 cat("CGEL of type: ", x$typeDesc, " (alpha = ", x$CGEL, ")\n") 192 193 if (!is.null(x$smooth)) 194 { 195 cat("Kernel: ", x$smooth$kernel," (bw=", x$smooth$bw,")\n\n") 196 }else { 197 cat("\n") 198 } 199 200 cat("Coefficients:\n") 201 print.default(format(x$coefficients, digits = digits), 202 print.gap = 2, quote = FALSE) 203 204 if (length(x$coefficients)<length(x$lambda)) 205 { 206 cat("\nLambdas:\n") 207 print.default(format(x$lambda, digits=digits), 208 print.gap = 2, quote = FALSE) 209 } else { 210 cat("\nNo table for Lambda when the model is just identified\n") 211 } 212 cat("\n", x$stest$ntest, "\n") 213 print.default(format(x$stest$test, digits=digits), 214 print.gap = 2, quote = FALSE) 215 cat("\n",x$specMod) 216 cat("\nConvergence code for the coefficients: ", x$conv_par, "\n") 217 cat("\nConvergence code for the lambdas: ", x$conv_lambda$convergence, "\n") 218 219 invisible(x) 220 } 221 222summary.gel <- function(object, ...) 223 { 224 z <- object 225 n <- nrow(z$gt) 226 se_par <- sqrt(diag(z$vcov_par)) 227 par <- z$coefficients 228 tval <- par/se_par 229 230 se_parl <- sqrt(diag(z$vcov_lambda)) 231 lamb <- z$lambda 232 tvall <- lamb/se_parl 233 234 ans <- list(type = z$type, call = z$call) 235 names(ans$type) <-"Type of GEL" 236 237 ans$coefficients <- round(cbind(par, se_par, tval, 2 * pnorm(abs(tval), lower.tail = FALSE)), 5) 238 ans$lambda <- round(cbind(lamb,se_parl, tvall, 2 * pnorm(abs(tvall), lower.tail = FALSE)), 5) 239 240 dimnames(ans$coefficients) <- list(names(z$coefficients), 241 c("Estimate", "Std. Error", "t value", "Pr(>|t|)")) 242 dimnames(ans$lambda) <- list(names(z$lambda), 243 c("Estimate", "Std. Error", "t value", "Pr(>|t|)")) 244 245 ans$stest=specTest(z) 246 247 if (z$type == "EL") 248 ans$badrho <- z$badrho 249 if (!is.null(z$weights)) 250 { 251 ans$weights <- z$weights 252 } 253 ans$conv_par <- z$conv_par 254 ans$conv_pt <- z$conv_pt 255 ans$conv_moment <- cbind(z$conv_moment) 256 ans$conv_lambda <- z$conv_lambda 257 ans$CGEL <- z$CGEL 258 ans$typeDesc <- z$typeDesc 259 ans$specMod <- z$specMod 260 if (!is.null(attr(object$dat,"smooth"))) 261 ans$smooth <- attr(object$dat,"smooth") 262 names(ans$conv_pt) <- "Sum_of_pt" 263 dimnames(ans$conv_moment) <- list(names(z$gt), "Sample_moment_with_pt") 264 class(ans) <- "summary.gel" 265 ans 266} 267 268residuals.gel <- function(object, ...) 269 { 270 if(is.null(object$model)) 271 stop("The residuals method is valid only for g=formula") 272 object$residuals 273 } 274 275fitted.gel <- function(object, ...) 276 { 277 if(is.null(object$model)) 278 stop("The residuals method is valid only for g=formula") 279 object$fitted.value 280 } 281 282formula.gel <- function(x, ...) 283{ 284 if(is.null(x$terms)) 285 stop("The gel object was not created by a formula") 286 else 287 formula(x$terms) 288} 289 290estfun.gel <- function(x, ...) 291 { 292 stop("estfun is not yet available for gel objects") 293 } 294 295bread.gel <- function(x, ...) 296 { 297 stop("Bread is not yet available for gel objects") 298 } 299 300 301getImpProb <- function(object, ...) 302 UseMethod("getImpProb") 303 304getImpProb.gel <- function(object, posProb=TRUE, normalize=TRUE, 305 checkConv=FALSE, ...) 306 { 307 if (!normalize || (object$type == "CUE" && !posProb)) 308 { 309 n <- NROW(object$gt) 310 pt <- -.rho(object$gt, object$lambda, type = object$type, 311 derive = 1, k = object$k1/object$k2)/n 312 if (object$type == "CUE" && posProb) 313 { 314 eps <- -length(pt)*min(min(pt),0) 315 pt <- (pt+eps/length(pt))/(1+eps) 316 } 317 if (normalize) 318 pt <- pt/sum(pt) 319 } else { 320 pt <- object$pt 321 } 322 if (checkConv) 323 attr(pt, "convergence") <- list(pt=sum(pt), 324 ptgt=colSums(pt*as.matrix(object$gt))) 325 pt 326 } 327 328 329 330