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