Nothing
R/SUMMARY.R
In gamlss: Generalised Additive Models for Location Scale and Shape
# this is a test summary for using the vcov matrix
# it should replace he old summary.gamlss
#---------------------------------------------------------------------------------------
summary.gamlss<- function (object,
type = c("vcov", "qr"), robust = FALSE,
save = FALSE,
hessian.fun = c("R", "PB"),
digits = max(3, getOption("digits") - 3),
...)
{
type <- match.arg(type)
## if use the vcov for se's
if (type=="vcov")
{
## this will be moved
covmat <- try(suppressWarnings(vcov(object, type="all", robust=robust, hessian.fun = hessian.fun)), silent = TRUE)
if (any(class(covmat)%in%"try-error"||any(is.na(covmat$se))))
{
warning(paste("summary: vcov has failed, option qr is used instead\n"))
type <- "qr"
}
## if try fails call old summary
}
ifWarning <- rep(FALSE, length(object$parameters))# to create warnings
if (type=="vcov")# type vcov --------------------------------------------------
{
coef <- covmat$coef
se <- covmat$se
tvalue <- coef/se
pvalue <- 2 * pt(-abs(tvalue), object$df.res) #if (est.disp) 2 * pt(-abs(tvalue), df.r) else 2 * pnorm(-abs(tvalue))
coef.table <- cbind(coef, se, tvalue, pvalue)
dimnames(coef.table) <- list(names(coef), c("Estimate" , "Std. Error" ,"t value","Pr(>|t|)"))
#now the table contains all informatio I need
#printCoefmat(coef.table, digits = digits, signif.stars = TRUE)
# to print use coef.table[1:3,]
cat("******************************************************************")
cat("\nFamily: ", deparse(object$family), "\n")
cat("\nCall: ", deparse(object$call,width.cutoff=50), "\n", fill=TRUE)
cat("Fitting method:", deparse(object$method), "\n\n")
est.disp <- FALSE
# df.r <- object$noObs - object$mu.df
#================ mu ESTIMATES ========================
if ("mu"%in%object$parameters)
{
# if (length(eval(parse(text=paste(paste("object$","mu", sep=""),".fix==TRUE", sep=""))))!=0)
if (length(eval(parse(text=paste(paste("object$","mu", sep=""),".fix==TRUE", sep=""))))!=0)
{
pm <- 1
p1 <- 1:pm
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Mu parameter is fixed \n")
if (all(object$mu.fv == object$mu.fv[1]))
cat("Mu = ", object$mu.fv[1], "\n")
else
cat("Mu is equal with the vector (", object$mu.fv[1], ",",object$mu.fv[2], ",",object$mu.fv[3], ",",object$mu.fv[4], ", ...) \n")
} else
{
ifWarning[1] <- (!is.null(unlist(attr(terms(formula(object, "mu"),
specials = .gamlss.sm.list), "specials"))))
if (object$mu.df != 0)
{
pm <- object$mu.qr$rank
p1 <- 1:pm
cat("------------------------------------------------------------------\n")
cat("Mu link function: ", object$mu.link)
cat("\n")
cat("Mu Coefficients:")
if (is.character(co <- object$contrasts))
cat(" [contrasts: ", apply(cbind(names(co), co), 1,
paste, collapse = "="), "]")
cat("\n")
printCoefmat(coef.table[p1,,drop=FALSE], digits = digits, signif.stars = TRUE)
# print.default(coef.table[p1,], digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
}
}
if ("sigma"%in%object$parameters)
{
if (length(eval(parse(text=paste(paste("object$","sigma", sep=""),".fix==TRUE", sep=""))))!=0)
{
ps <- 1
p1 <- (pm+1):(pm+ps)
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Sigma parameter is fixed \n")
if (all(object$sigma.fv == object$sigma.fv[1]))
cat("Sigma = ", object$sigma.fv[1], "\n")
else
cat("Sigma is equal with the vector (", object$sigma.fv[1], ",",object$sigma.fv[2], ",",object$sigma.fv[3], ",",object$sigma.fv[4], ", ...) \n")
} else
{
ifWarning[2] <- (!is.null(unlist(attr(terms(formula(object, "sigma"),
specials = .gamlss.sm.list), "specials"))))
if (object$sigma.df != 0)
{
ps <- object$sigma.qr$rank
p1 <- (pm+1):(pm+ps)
#cat("\n")
cat("------------------------------------------------------------------\n")
cat("Sigma link function: ", object$sigma.link)
cat("\n")
cat("Sigma Coefficients:")
cat("\n")
printCoefmat(coef.table[p1,, drop=FALSE], digits = digits, signif.stars = TRUE)
# print.default(coef.table[p1,], digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
}
}
if ("nu"%in%object$parameters)
{
if (length(eval(parse(text=paste(paste("object$","nu", sep=""),".fix==TRUE", sep=""))))!=0)
{
pn <- 1
p1 <- (pm+ps+1):(pm+ps+pn)
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Nu parameter is fixed \n")
if (all(object$nu.fv == object$nu.fv[1]))
cat("Nu = ", object$nu.fv[1], "\n")
else
cat("nu is equal with the vector (", object$nu.fv[1], ",",object$nu.fv[2], ",",object$nu.fv[3], ",",object$nu.fv[4], ", ...) \n")
} else
{
ifWarning[3] <- (!is.null(unlist(attr(terms(formula(object, "nu"),
specials = .gamlss.sm.list), "specials"))))
if (object$nu.df != 0)
{
pn <- object$nu.qr$rank
p1 <- (pm+ps+1):(pm+ps+pn)
cat("------------------------------------------------------------------\n")
cat("Nu link function: ", object$nu.link,"\n")
cat("Nu Coefficients:")
cat("\n")
printCoefmat(coef.table[p1,, drop=FALSE], digits = digits, signif.stars = TRUE)
# print.default(coef.table[p1,], digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
}
}
if ("tau"%in%object$parameters)
{
if (length(eval(parse(text=paste(paste("object$","tau", sep=""),".fix==TRUE", sep=""))))!=0)
{
pt <- 1
p1 <- (pm+ps+pn+1):(pm+ps+pn+pt)
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Tau parameter is fixed \n")
if (all(object$tau.fv == object$tau.fv[1]))
cat("Tau = ", object$tau.fv[1], "\n")
else
cat("Tau is equal with the vector (", object$tau.fv[1], ",",object$tau.fv[2], ",",object$tau.fv[3], ",",object$tau.fv[4], ", ...) \n")
} else
{
ifWarning[4] <- (!is.null(unlist(attr(terms(formula(object, "tau"),
specials = .gamlss.sm.list), "specials"))))
if (object$tau.df != 0)
{
pt <- object$tau.qr$rank
p1 <- (pm+ps+pn+1):(pm+ps+pn+pt)
cat("------------------------------------------------------------------\n")
cat("Tau link function: ", object$tau.link,"\n")
cat("Tau Coefficients:")
cat("\n")
printCoefmat(coef.table[p1,, drop=FALSE], digits = digits, signif.stars = TRUE)
# print.default(coef.table[p1,], digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
}
}
if (any(ifWarning))
{
cat("------------------------------------------------------------------\n")
cat("NOTE: Additive smoothing terms exist in the formulas: \n")
cat(" i) Std. Error for smoothers are for the linear effect only. \n")
cat("ii) Std. Error for the linear terms maybe are not accurate. \n")
}
cat("------------------------------------------------------------------\n")
cat("No. of observations in the fit: ", object$noObs, "\n")
cat("Degrees of Freedom for the fit: ", object$df.fit)
cat("\n")
cat(" Residual Deg. of Freedom: ", object$df.residual, "\n")
cat(" at cycle: ", object$iter, "\n \n")
cat("Global Deviance: ", object$G.deviance,#format(signif(object$G.deviance, digits)),
"\n AIC: ",object$aic, #format(signif(object$aic, digits)),
"\n SBC: ",object$sbc, "\n") #format(signif(object$sbc, digits)), "\n")
cat("******************************************************************")
cat("\n")
}
if (type=="qr")# TYPE qr ---------------------------------------------------
{
# local function definition
#-------------------------------------------------------------------------------
estimatesgamlss<-function (object,Qr, p1, coef.p, est.disp , df.r, digits = max(3, getOption("digits") - 3),
covmat.unscaled , ...)
{
#covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
dimnames(covmat.unscaled) <- list(names(coef.p), names(coef.p))
covmat <- covmat.unscaled #in glm is=dispersion * covmat.unscaled, but here is already multiplied by the dispersion
var.cf <- diag(covmat)
s.err <- sqrt(var.cf)
tvalue <- coef.p/s.err
dn <- c("Estimate", "Std. Error")
if (!est.disp)
{
pvalue <- 2 * pnorm(-abs(tvalue))
coef.table <- cbind(coef.p, s.err, tvalue, pvalue)
dimnames(coef.table) <- list(names(coef.p), c(dn, "z value",
"Pr(>|z|)"))
}
else if (df.r > 0)
{
pvalue <- 2 * pt(-abs(tvalue), df.r)
coef.table <- cbind(coef.p, s.err, tvalue, pvalue)
dimnames(coef.table) <- list(names(coef.p), c(dn, "t value","Pr(>|t|)"))
}
else
{
coef.table <- cbind(coef.p, Inf)
dimnames(coef.table) <- list(names(coef.p), dn)
}
return(coef.table)
}
##------------------------------------------------------------------------------
# here for the proper function
##------------------------------------------------------------------------------
dispersion <- NULL
# digits <- max(3, getOption("digits") - 3)
cat("******************************************************************")
cat("\nFamily: ", deparse(object$family), "\n")
cat("\nCall: ", deparse(object$call), "\n", fill=TRUE)
cat("Fitting method:", deparse(object$method), "\n\n")
est.disp <- FALSE
df.r <- object$noObs - object$mu.df
#================ mu ESTIMATES ========================
if ("mu"%in%object$parameters)
{
ifWarning[1] <- (!is.null(unlist(attr(terms(formula(object, "mu"),
specials = .gamlss.sm.list), "specials"))))
if (object$mu.df != 0)
{
Qr <- object$mu.qr
df.r <- object$noObs - object$mu.df
## this should be taken out sinse there is no dispersion MS Thursday, August 17, 2006
if (is.null(dispersion))
dispersion <- if (any(object$family == c("PO",
"BI", "EX", "P1")))
1
else if (df.r > 0) {
est.disp <- TRUE
if (any(object$weights == 0))
warning(paste("observations with zero weight",
"not used for calculating dispersion"))
}
else Inf
#---------------------------------------------------end of taken out
p <- object$mu.df # object$rank
p1 <- 1:(p-object$mu.nl.df)
covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
mu.coef.table <- estimatesgamlss(object=object,Qr=object$mu.qr, p1=p1,
coef.p=object$mu.coefficients[Qr$pivot[p1]],
est.disp =est.disp, df.r=df.r,
covmat.unscaled =covmat.unscaled )
cat("------------------------------------------------------------------\n")
cat("Mu link function: ", object$mu.link)
cat("\n")
cat("Mu Coefficients:")
if (is.character(co <- object$contrasts))
cat(" [contrasts: ", apply(cbind(names(co), co), 1,
paste, collapse = "="), "]")
cat("\n")
printCoefmat(mu.coef.table, digits = digits, signif.stars = TRUE)
# print.default(mu.coef.table, digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
else
if(object$mu.fix == TRUE)
{
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Mu parameter is fixed")
cat("\n")
if (all(object$mu.fv == object$mu.fv[1]))
cat("Mu = ", object$mu.fv[1], "\n")
else
cat("Mu is equal with the vector (", object$mu.fv[1], ",",object$mu.fv[2], ",",object$mu.fv[3], ",",object$mu.fv[4], ", ...) \n")
}
coef.table <- mu.coef.table
}
else
{
if (df.r > 0) {
est.disp <- TRUE
if (any(object$weights == 0))
warning(paste("observations with zero weight",
"not used for calculating dispersion"))
}
}
if ("sigma"%in%object$parameters)
{
ifWarning[2] <- (!is.null(unlist(attr(terms(formula(object, "sigma"),
specials = .gamlss.sm.list), "specials"))))
if (object$sigma.df != 0)
{
Qr <- object$sigma.qr
df.r <- object$noObs - object$sigma.df
p <- object$sigma.df
p1 <- 1:(p-object$sigma.nl.df)
covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
sigma.coef.table<-estimatesgamlss(object=object,Qr=object$sigma.qr, p1=p1,
coef.p=object$sigma.coefficients[Qr$pivot[p1]],
est.disp =est.disp, df.r=df.r,
covmat.unscaled =covmat.unscaled )
#cat("\n")
cat("------------------------------------------------------------------\n")
cat("Sigma link function: ", object$sigma.link)
cat("\n")
cat("Sigma Coefficients:")
cat("\n")
printCoefmat(sigma.coef.table, digits = digits, signif.stars = TRUE)
#print.default(sigma.coef.table, digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
else
if(object$sigma.fix == TRUE) {
cat("------------------------------------------------------------------\n")
cat("Sigma parameter is fixed")
cat("\n")
if (all(object$sigma.fv == object$sigma.fv[1]))
cat("Sigma = ", object$sigma.fv[1], "\n")
else
cat("Sigma is equal with the vector (", object$sigma.fv[1], ",",object$sigma.fv[2], ",",object$sigma.fv[3], ",",object$sigma.fv[4], ", ...) \n")
}
coef.table <- rbind(mu.coef.table, sigma.coef.table)
}
if ("nu"%in%object$parameters)
{
ifWarning[3] <- (!is.null(unlist(attr(terms(formula(object, "nu"),
specials = .gamlss.sm.list), "specials"))))
if (object$nu.df != 0)
{
Qr <- object$nu.qr
df.r <- object$noObs - object$nu.df
p <- object$nu.df
p1 <- 1:(p-object$nu.nl.df)
covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
nu.coef.table <- estimatesgamlss(object=object,Qr=object$nu.qr, p1=p1,
coef.p=object$nu.coefficients[Qr$pivot[p1]],
est.disp =est.disp, df.r=df.r,
covmat.unscaled =covmat.unscaled )
cat("------------------------------------------------------------------\n")
cat("Nu link function: ", object$nu.link,"\n")
cat("Nu Coefficients:")
cat("\n")
printCoefmat(nu.coef.table, digits = digits, signif.stars = TRUE)
#print.default(nu.coef.table, digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
else
if(object$nu.fix == TRUE) {
cat("------------------------------------------------------------------\n")
cat("Nu parameter is fixed")
cat("\n")
if (all(object$nu.fv == object$nu.fv[1]))
cat("Nu = ", object$nu.fv[1], "\n")
else
cat("Nu is equal with the vector (", object$nu.fv[1], ",",object$nu.fv[2], ",",object$nu.fv[3], ",",object$nu.fv[4], ", ...) \n")
}
coef.table <- rbind(mu.coef.table, sigma.coef.table, nu.coef.table)
}
if ("tau"%in%object$parameters)
{
ifWarning[4] <- (!is.null(unlist(attr(terms(formula(object, "tau"),
specials = .gamlss.sm.list), "specials"))))
if (object$tau.df != 0)
{
Qr <- object$tau.qr
df.r <- object$noObs - object$tau.df
p <- object$tau.df
p1 <- 1:(p-object$tau.nl.df)
covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
tau.coef.table <- estimatesgamlss(object=object,Qr=object$tau.qr, p1=p1,
coef.p=object$tau.coefficients[Qr$pivot[p1]],
est.disp =est.disp, df.r=df.r,
covmat.unscaled =covmat.unscaled )
cat("------------------------------------------------------------------\n")
cat("Tau link function: ", object$tau.link,"\n")
cat("Tau Coefficients:")
cat("\n")
printCoefmat(tau.coef.table, digits = digits, signif.stars = TRUE)
#print.default(tau.coef.table, digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
else
if(object$tau.fix == TRUE) {
cat("------------------------------------------------------------------\n")
cat("Tau parameter is fixed")
cat("\n")
if (all(object$tau.fv == object$tau.fv[1]))
cat("Tau = ", object$tau.fv[1], "\n")
else
cat("Tau is equal with the vector (", object$tau.fv[1], ",",object$tau.fv[2], ",",object$tau.fv[3], ",",object$tau.fv[4], ", ...) \n")
}
coef.table <- rbind(mu.coef.table, sigma.coef.table, nu.coef.table, tau.coef.table)
}
if (any(ifWarning))
{
cat("------------------------------------------------------------------\n")
cat("NOTE: Additive smoothing terms exist in the formulas: \n")
cat(" i) Std. Error for smoothers are for the linear effect only. \n")
cat("ii) Std. Error for the linear terms may not be reliable. \n")
}
cat("------------------------------------------------------------------\n")
cat("No. of observations in the fit: ", object$noObs, "\n")
cat("Degrees of Freedom for the fit: ", object$df.fit)
cat("\n")
cat(" Residual Deg. of Freedom: ", object$df.residual, "\n")
cat(" at cycle: ", object$iter, "\n \n")
cat("Global Deviance: ", object$G.deviance,#format(signif(object$G.deviance, digits)),
"\n AIC: ",object$aic, #format(signif(object$aic, digits)),
"\n SBC: ",object$sbc, "\n") #format(signif(object$sbc, digits)), "\n")
cat("******************************************************************")
cat("\n")
}
if ( save == TRUE)
{
out <- as.list(environment())
return(out)
}
invisible(coef.table)
}
#==============================================================================
Try the gamlss package in your browser
Any scripts or data that you put into this service are public.
gamlss documentation built on Oct. 4, 2023, 5:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# this is a test summary for using the vcov matrix
# it should replace he old summary.gamlss
#---------------------------------------------------------------------------------------
summary.gamlss<- function (object,
type = c("vcov", "qr"), robust = FALSE,
save = FALSE,
hessian.fun = c("R", "PB"),
digits = max(3, getOption("digits") - 3),
...)
{
type <- match.arg(type)
## if use the vcov for se's
if (type=="vcov")
{
## this will be moved
covmat <- try(suppressWarnings(vcov(object, type="all", robust=robust, hessian.fun = hessian.fun)), silent = TRUE)
if (any(class(covmat)%in%"try-error"||any(is.na(covmat$se))))
{
warning(paste("summary: vcov has failed, option qr is used instead\n"))
type <- "qr"
}
## if try fails call old summary
}
ifWarning <- rep(FALSE, length(object$parameters))# to create warnings
if (type=="vcov")# type vcov --------------------------------------------------
{
coef <- covmat$coef
se <- covmat$se
tvalue <- coef/se
pvalue <- 2 * pt(-abs(tvalue), object$df.res) #if (est.disp) 2 * pt(-abs(tvalue), df.r) else 2 * pnorm(-abs(tvalue))
coef.table <- cbind(coef, se, tvalue, pvalue)
dimnames(coef.table) <- list(names(coef), c("Estimate" , "Std. Error" ,"t value","Pr(>|t|)"))
#now the table contains all informatio I need
#printCoefmat(coef.table, digits = digits, signif.stars = TRUE)
# to print use coef.table[1:3,]
cat("******************************************************************")
cat("\nFamily: ", deparse(object$family), "\n")
cat("\nCall: ", deparse(object$call,width.cutoff=50), "\n", fill=TRUE)
cat("Fitting method:", deparse(object$method), "\n\n")
est.disp <- FALSE
# df.r <- object$noObs - object$mu.df
#================ mu ESTIMATES ========================
if ("mu"%in%object$parameters)
{
# if (length(eval(parse(text=paste(paste("object$","mu", sep=""),".fix==TRUE", sep=""))))!=0)
if (length(eval(parse(text=paste(paste("object$","mu", sep=""),".fix==TRUE", sep=""))))!=0)
{
pm <- 1
p1 <- 1:pm
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Mu parameter is fixed \n")
if (all(object$mu.fv == object$mu.fv[1]))
cat("Mu = ", object$mu.fv[1], "\n")
else
cat("Mu is equal with the vector (", object$mu.fv[1], ",",object$mu.fv[2], ",",object$mu.fv[3], ",",object$mu.fv[4], ", ...) \n")
} else
{
ifWarning[1] <- (!is.null(unlist(attr(terms(formula(object, "mu"),
specials = .gamlss.sm.list), "specials"))))
if (object$mu.df != 0)
{
pm <- object$mu.qr$rank
p1 <- 1:pm
cat("------------------------------------------------------------------\n")
cat("Mu link function: ", object$mu.link)
cat("\n")
cat("Mu Coefficients:")
if (is.character(co <- object$contrasts))
cat(" [contrasts: ", apply(cbind(names(co), co), 1,
paste, collapse = "="), "]")
cat("\n")
printCoefmat(coef.table[p1,,drop=FALSE], digits = digits, signif.stars = TRUE)
# print.default(coef.table[p1,], digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
}
}
if ("sigma"%in%object$parameters)
{
if (length(eval(parse(text=paste(paste("object$","sigma", sep=""),".fix==TRUE", sep=""))))!=0)
{
ps <- 1
p1 <- (pm+1):(pm+ps)
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Sigma parameter is fixed \n")
if (all(object$sigma.fv == object$sigma.fv[1]))
cat("Sigma = ", object$sigma.fv[1], "\n")
else
cat("Sigma is equal with the vector (", object$sigma.fv[1], ",",object$sigma.fv[2], ",",object$sigma.fv[3], ",",object$sigma.fv[4], ", ...) \n")
} else
{
ifWarning[2] <- (!is.null(unlist(attr(terms(formula(object, "sigma"),
specials = .gamlss.sm.list), "specials"))))
if (object$sigma.df != 0)
{
ps <- object$sigma.qr$rank
p1 <- (pm+1):(pm+ps)
#cat("\n")
cat("------------------------------------------------------------------\n")
cat("Sigma link function: ", object$sigma.link)
cat("\n")
cat("Sigma Coefficients:")
cat("\n")
printCoefmat(coef.table[p1,, drop=FALSE], digits = digits, signif.stars = TRUE)
# print.default(coef.table[p1,], digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
}
}
if ("nu"%in%object$parameters)
{
if (length(eval(parse(text=paste(paste("object$","nu", sep=""),".fix==TRUE", sep=""))))!=0)
{
pn <- 1
p1 <- (pm+ps+1):(pm+ps+pn)
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Nu parameter is fixed \n")
if (all(object$nu.fv == object$nu.fv[1]))
cat("Nu = ", object$nu.fv[1], "\n")
else
cat("nu is equal with the vector (", object$nu.fv[1], ",",object$nu.fv[2], ",",object$nu.fv[3], ",",object$nu.fv[4], ", ...) \n")
} else
{
ifWarning[3] <- (!is.null(unlist(attr(terms(formula(object, "nu"),
specials = .gamlss.sm.list), "specials"))))
if (object$nu.df != 0)
{
pn <- object$nu.qr$rank
p1 <- (pm+ps+1):(pm+ps+pn)
cat("------------------------------------------------------------------\n")
cat("Nu link function: ", object$nu.link,"\n")
cat("Nu Coefficients:")
cat("\n")
printCoefmat(coef.table[p1,, drop=FALSE], digits = digits, signif.stars = TRUE)
# print.default(coef.table[p1,], digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
}
}
if ("tau"%in%object$parameters)
{
if (length(eval(parse(text=paste(paste("object$","tau", sep=""),".fix==TRUE", sep=""))))!=0)
{
pt <- 1
p1 <- (pm+ps+pn+1):(pm+ps+pn+pt)
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Tau parameter is fixed \n")
if (all(object$tau.fv == object$tau.fv[1]))
cat("Tau = ", object$tau.fv[1], "\n")
else
cat("Tau is equal with the vector (", object$tau.fv[1], ",",object$tau.fv[2], ",",object$tau.fv[3], ",",object$tau.fv[4], ", ...) \n")
} else
{
ifWarning[4] <- (!is.null(unlist(attr(terms(formula(object, "tau"),
specials = .gamlss.sm.list), "specials"))))
if (object$tau.df != 0)
{
pt <- object$tau.qr$rank
p1 <- (pm+ps+pn+1):(pm+ps+pn+pt)
cat("------------------------------------------------------------------\n")
cat("Tau link function: ", object$tau.link,"\n")
cat("Tau Coefficients:")
cat("\n")
printCoefmat(coef.table[p1,, drop=FALSE], digits = digits, signif.stars = TRUE)
# print.default(coef.table[p1,], digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
}
}
if (any(ifWarning))
{
cat("------------------------------------------------------------------\n")
cat("NOTE: Additive smoothing terms exist in the formulas: \n")
cat(" i) Std. Error for smoothers are for the linear effect only. \n")
cat("ii) Std. Error for the linear terms maybe are not accurate. \n")
}
cat("------------------------------------------------------------------\n")
cat("No. of observations in the fit: ", object$noObs, "\n")
cat("Degrees of Freedom for the fit: ", object$df.fit)
cat("\n")
cat(" Residual Deg. of Freedom: ", object$df.residual, "\n")
cat(" at cycle: ", object$iter, "\n \n")
cat("Global Deviance: ", object$G.deviance,#format(signif(object$G.deviance, digits)),
"\n AIC: ",object$aic, #format(signif(object$aic, digits)),
"\n SBC: ",object$sbc, "\n") #format(signif(object$sbc, digits)), "\n")
cat("******************************************************************")
cat("\n")
}
if (type=="qr")# TYPE qr ---------------------------------------------------
{
# local function definition
#-------------------------------------------------------------------------------
estimatesgamlss<-function (object,Qr, p1, coef.p, est.disp , df.r, digits = max(3, getOption("digits") - 3),
covmat.unscaled , ...)
{
#covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
dimnames(covmat.unscaled) <- list(names(coef.p), names(coef.p))
covmat <- covmat.unscaled #in glm is=dispersion * covmat.unscaled, but here is already multiplied by the dispersion
var.cf <- diag(covmat)
s.err <- sqrt(var.cf)
tvalue <- coef.p/s.err
dn <- c("Estimate", "Std. Error")
if (!est.disp)
{
pvalue <- 2 * pnorm(-abs(tvalue))
coef.table <- cbind(coef.p, s.err, tvalue, pvalue)
dimnames(coef.table) <- list(names(coef.p), c(dn, "z value",
"Pr(>|z|)"))
}
else if (df.r > 0)
{
pvalue <- 2 * pt(-abs(tvalue), df.r)
coef.table <- cbind(coef.p, s.err, tvalue, pvalue)
dimnames(coef.table) <- list(names(coef.p), c(dn, "t value","Pr(>|t|)"))
}
else
{
coef.table <- cbind(coef.p, Inf)
dimnames(coef.table) <- list(names(coef.p), dn)
}
return(coef.table)
}
##------------------------------------------------------------------------------
# here for the proper function
##------------------------------------------------------------------------------
dispersion <- NULL
# digits <- max(3, getOption("digits") - 3)
cat("******************************************************************")
cat("\nFamily: ", deparse(object$family), "\n")
cat("\nCall: ", deparse(object$call), "\n", fill=TRUE)
cat("Fitting method:", deparse(object$method), "\n\n")
est.disp <- FALSE
df.r <- object$noObs - object$mu.df
#================ mu ESTIMATES ========================
if ("mu"%in%object$parameters)
{
ifWarning[1] <- (!is.null(unlist(attr(terms(formula(object, "mu"),
specials = .gamlss.sm.list), "specials"))))
if (object$mu.df != 0)
{
Qr <- object$mu.qr
df.r <- object$noObs - object$mu.df
## this should be taken out sinse there is no dispersion MS Thursday, August 17, 2006
if (is.null(dispersion))
dispersion <- if (any(object$family == c("PO",
"BI", "EX", "P1")))
1
else if (df.r > 0) {
est.disp <- TRUE
if (any(object$weights == 0))
warning(paste("observations with zero weight",
"not used for calculating dispersion"))
}
else Inf
#---------------------------------------------------end of taken out
p <- object$mu.df # object$rank
p1 <- 1:(p-object$mu.nl.df)
covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
mu.coef.table <- estimatesgamlss(object=object,Qr=object$mu.qr, p1=p1,
coef.p=object$mu.coefficients[Qr$pivot[p1]],
est.disp =est.disp, df.r=df.r,
covmat.unscaled =covmat.unscaled )
cat("------------------------------------------------------------------\n")
cat("Mu link function: ", object$mu.link)
cat("\n")
cat("Mu Coefficients:")
if (is.character(co <- object$contrasts))
cat(" [contrasts: ", apply(cbind(names(co), co), 1,
paste, collapse = "="), "]")
cat("\n")
printCoefmat(mu.coef.table, digits = digits, signif.stars = TRUE)
# print.default(mu.coef.table, digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
else
if(object$mu.fix == TRUE)
{
#warning(paste("mu parameter is fixed"))
cat("------------------------------------------------------------------\n")
cat("Mu parameter is fixed")
cat("\n")
if (all(object$mu.fv == object$mu.fv[1]))
cat("Mu = ", object$mu.fv[1], "\n")
else
cat("Mu is equal with the vector (", object$mu.fv[1], ",",object$mu.fv[2], ",",object$mu.fv[3], ",",object$mu.fv[4], ", ...) \n")
}
coef.table <- mu.coef.table
}
else
{
if (df.r > 0) {
est.disp <- TRUE
if (any(object$weights == 0))
warning(paste("observations with zero weight",
"not used for calculating dispersion"))
}
}
if ("sigma"%in%object$parameters)
{
ifWarning[2] <- (!is.null(unlist(attr(terms(formula(object, "sigma"),
specials = .gamlss.sm.list), "specials"))))
if (object$sigma.df != 0)
{
Qr <- object$sigma.qr
df.r <- object$noObs - object$sigma.df
p <- object$sigma.df
p1 <- 1:(p-object$sigma.nl.df)
covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
sigma.coef.table<-estimatesgamlss(object=object,Qr=object$sigma.qr, p1=p1,
coef.p=object$sigma.coefficients[Qr$pivot[p1]],
est.disp =est.disp, df.r=df.r,
covmat.unscaled =covmat.unscaled )
#cat("\n")
cat("------------------------------------------------------------------\n")
cat("Sigma link function: ", object$sigma.link)
cat("\n")
cat("Sigma Coefficients:")
cat("\n")
printCoefmat(sigma.coef.table, digits = digits, signif.stars = TRUE)
#print.default(sigma.coef.table, digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
else
if(object$sigma.fix == TRUE) {
cat("------------------------------------------------------------------\n")
cat("Sigma parameter is fixed")
cat("\n")
if (all(object$sigma.fv == object$sigma.fv[1]))
cat("Sigma = ", object$sigma.fv[1], "\n")
else
cat("Sigma is equal with the vector (", object$sigma.fv[1], ",",object$sigma.fv[2], ",",object$sigma.fv[3], ",",object$sigma.fv[4], ", ...) \n")
}
coef.table <- rbind(mu.coef.table, sigma.coef.table)
}
if ("nu"%in%object$parameters)
{
ifWarning[3] <- (!is.null(unlist(attr(terms(formula(object, "nu"),
specials = .gamlss.sm.list), "specials"))))
if (object$nu.df != 0)
{
Qr <- object$nu.qr
df.r <- object$noObs - object$nu.df
p <- object$nu.df
p1 <- 1:(p-object$nu.nl.df)
covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
nu.coef.table <- estimatesgamlss(object=object,Qr=object$nu.qr, p1=p1,
coef.p=object$nu.coefficients[Qr$pivot[p1]],
est.disp =est.disp, df.r=df.r,
covmat.unscaled =covmat.unscaled )
cat("------------------------------------------------------------------\n")
cat("Nu link function: ", object$nu.link,"\n")
cat("Nu Coefficients:")
cat("\n")
printCoefmat(nu.coef.table, digits = digits, signif.stars = TRUE)
#print.default(nu.coef.table, digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
else
if(object$nu.fix == TRUE) {
cat("------------------------------------------------------------------\n")
cat("Nu parameter is fixed")
cat("\n")
if (all(object$nu.fv == object$nu.fv[1]))
cat("Nu = ", object$nu.fv[1], "\n")
else
cat("Nu is equal with the vector (", object$nu.fv[1], ",",object$nu.fv[2], ",",object$nu.fv[3], ",",object$nu.fv[4], ", ...) \n")
}
coef.table <- rbind(mu.coef.table, sigma.coef.table, nu.coef.table)
}
if ("tau"%in%object$parameters)
{
ifWarning[4] <- (!is.null(unlist(attr(terms(formula(object, "tau"),
specials = .gamlss.sm.list), "specials"))))
if (object$tau.df != 0)
{
Qr <- object$tau.qr
df.r <- object$noObs - object$tau.df
p <- object$tau.df
p1 <- 1:(p-object$tau.nl.df)
covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
tau.coef.table <- estimatesgamlss(object=object,Qr=object$tau.qr, p1=p1,
coef.p=object$tau.coefficients[Qr$pivot[p1]],
est.disp =est.disp, df.r=df.r,
covmat.unscaled =covmat.unscaled )
cat("------------------------------------------------------------------\n")
cat("Tau link function: ", object$tau.link,"\n")
cat("Tau Coefficients:")
cat("\n")
printCoefmat(tau.coef.table, digits = digits, signif.stars = TRUE)
#print.default(tau.coef.table, digits = digits, print.gap = 2, quote = FALSE)
cat("\n")
}
else
if(object$tau.fix == TRUE) {
cat("------------------------------------------------------------------\n")
cat("Tau parameter is fixed")
cat("\n")
if (all(object$tau.fv == object$tau.fv[1]))
cat("Tau = ", object$tau.fv[1], "\n")
else
cat("Tau is equal with the vector (", object$tau.fv[1], ",",object$tau.fv[2], ",",object$tau.fv[3], ",",object$tau.fv[4], ", ...) \n")
}
coef.table <- rbind(mu.coef.table, sigma.coef.table, nu.coef.table, tau.coef.table)
}
if (any(ifWarning))
{
cat("------------------------------------------------------------------\n")
cat("NOTE: Additive smoothing terms exist in the formulas: \n")
cat(" i) Std. Error for smoothers are for the linear effect only. \n")
cat("ii) Std. Error for the linear terms may not be reliable. \n")
}
cat("------------------------------------------------------------------\n")
cat("No. of observations in the fit: ", object$noObs, "\n")
cat("Degrees of Freedom for the fit: ", object$df.fit)
cat("\n")
cat(" Residual Deg. of Freedom: ", object$df.residual, "\n")
cat(" at cycle: ", object$iter, "\n \n")
cat("Global Deviance: ", object$G.deviance,#format(signif(object$G.deviance, digits)),
"\n AIC: ",object$aic, #format(signif(object$aic, digits)),
"\n SBC: ",object$sbc, "\n") #format(signif(object$sbc, digits)), "\n")
cat("******************************************************************")
cat("\n")
}
if ( save == TRUE)
{
out <- as.list(environment())
return(out)
}
invisible(coef.table)
}
#==============================================================================
Try the gamlss package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.