Nothing
R/gamlssMX-10-11-07.R
In gamlss.mx: Fitting Mixture Distributions with GAMLSS
#-------------------------------------------------------------------------------
## this a general function for fitting mixures within gamlss
## created by Mikis Stasinopoulos and Bob Rigby
## LAST CHECKED 5-8-14
## latest change Tuesday, April 10, 2007 at 08:54
## Some of the models which can be fitted
## with the function gamlssNP can be identical
## In this version the prior probabilities pi can be modelled
## TO DO
## it needs a summary functions but the EM-algorithm SE's are not correct
## do we need EBP here?
## do we need means?
## the plotMX() need modification for BI type data Done but not checked
##
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
gamlssMX <- function ( formula = formula(data),
pi.formula = ~1,
family = "NO", # note is character
weights,
K = 2,
prob = NULL,
data = sys.parent(),
control = MX.control(...),
g.control = gamlss.control(trace=FALSE,...),
zero.component = FALSE,
...)
{
#-------------------------------------------------------------------------------
# this is to replicate rqres within gamlss enviroment DS Friday, March 31, 2006 at 10:30
rqres <- function (pfun = "pNO",
type = c("Continuous", "Discrete", "Mixed"),
censored = NULL,
ymin = NULL,
mass.p = NULL,
prob.mp = NULL,
y = y,
... )
{ }
body(rqres) <- eval(quote(body(rqres)), envir = getNamespace("gamlss"))
##------------------------------------------------------------------------------
##------------------------------------------------------------------------------
.gamlss.bi.list <- eval(quote(.gamlss.bi.list), envir = getNamespace("gamlss"))
# extra functions within -------------------------------------------------------
# for getting the commulative function------------------------------------------
get.the.p.function <- function(object, ...)
{
#gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
if (!is.gamlss(object)) stop(paste("This is not an gamlss object", "\n", ""))
fname <- object$family[[1]]
DistPar <- object$parameters
nopar <- length(DistPar)
dfun <- paste("p",fname,sep="")
# binomial denominators
switch(nopar,
{pfun <- if (fname%in%.gamlss.bi.list) eval(call(dfun,q=object$y, bd=object$bd, mu=fitted(object,"mu")))
else eval(call(dfun,q=object$y, mu=fitted(object,"mu")))},
{pfun <- if (fname%in%.gamlss.bi.list) eval(call(dfun,q=object$y, bd=object$bd, mu=fitted(object,"mu"),sigma=fitted(object,"sigma")))
else eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"))) },
{pfun <- if (fname%in%.gamlss.bi.list) eval(call(dfun,q=object$y, bd=object$bd, mu=fitted(object,"mu"),sigma=fitted(object,"sigma"), nu=fitted(object,"nu")))
else eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"), nu=fitted(object,"nu")))},
{pfun <- if (fname%in%.gamlss.bi.list) eval(call(dfun,q=object$y, bd=object$bd, mu=fitted(object,"mu"),sigma=fitted(object,"sigma"), nu=fitted(object,"nu"), tau=fitted(object,"tau")))
else eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"), nu=fitted(object,"nu"), tau=fitted(object,"tau")))})
pfun
}
#-------------------------------------------------------------------------------
get.likelihood <- function(obj)
{
#gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
if (!is.gamlss(obj)) stop(paste("This is not an gamlss object", "\n", ""))
fname <- obj$family[[1]]
DistPar <- obj$parameters
nopar <- length(DistPar)
dfun <- paste("d",fname,sep="")
switch(nopar,
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun,x=obj$y, bd=obj$bd, mu=fitted(obj,"mu")))
else eval(call(dfun, x=obj$y, mu=fitted(obj)))},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun,x=obj$y, bd=obj$bd, mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma") ))
else eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"))) },
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun,x=obj$y, bd=obj$bd, mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu") ))
else eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu")))},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun,x=obj$y, bd=obj$bd, mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu"), tau=fitted(obj,"tau") ))
else eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu"), tau=fitted(obj,"tau")))})
lik
}
#-------------------------------------------------------------------------------
expand.vc <- function(x,ni)
{
if (length(ni)==1)
{
if (ni<1) stop("ni should be greater than 1")
xx <- x
if (ni==1) return(xx)
for ( i in 2:ni) xx <- rbind(x,xx)
xx
}
else
{
n <- dim(x)[[1]]
c <- dim(x)[[2]]
xx <- NULL
for ( i in seq(1,n))
{
xx <- rbind(xx,matrix(rep(x[i,],ni[i]),ncol=c,byrow=TRUE))
}
xx
}
}
#-------------------------------------------------------------------------------
# the proper function starts here
#-------------------------------------------------------------------------------
#library(gamlss)
#library(nnet)
gamlssMXcall <- match.call() # the function call
## checking for NA in the data
if(!missing(data) & any(is.na(data)))
stop("The data contains NA's, use data = na.omit(mydata)")
## prob are the starting values of the pi's
if(zero.component == TRUE)
{ KK <- K+1
if (is.null(prob))
{
prob <- rep(1/KK, KK)
}
else
{
if (length(prob)!=KK) stop("the length of prob must be K+1(for zero component) ")
}
}
else
{
if (is.null(prob))
{
prob <- rep(1/K, K)
}
else
{
if (length(prob)!=K) stop("the length of prob must be K")
}
}
## expand the formula if not list
if (!is.list(formula))
{
allFormula <- vector("list",K)
for (i in 1:K) allFormula[[i]] <- formula
}
else allFormula <- formula
## expand the family if not list
if (!is.list(family))
{
allFamily <- vector("list",K)
for (i in 1:K) allFamily[[i]] <- family
}
else allFamily <- family
## check that the length for both the family and formula is K
if (length(allFamily)!=K) stop("the length of the list for the family is not equal the number of componets in the mixture")
if (length(allFormula)!=K) stop("the length of the list for the formula is not equal the number of componets in the mixture")
## get response for his length
Y <- model.extract(model.frame(allFormula[[1]],data=data),"response")
N <- if(is.null(dim(Y))) length(Y) else dim(Y)[1]# calculate the dimension for y
pweights. <- pweights. <<- if (missing(weights)) rep(1,N) else weights
## get things from control
set.seed(control$seed) # set the seed
prob.sample <- if (is.null(control$sample)) 10/N else control$sample
## creating the matrix of posterior probabilities (weights)
W <- matrix(1, nrow=N, ncol=K )
allModels <- vector("list",K)
## get starting values for the weights (posterior probabilities)
for (i in 1:K)
{ wSam <- sample(c(0,1), N, replace = TRUE, prob=c(1-min(.5,prob.sample),min(.5,prob.sample)))
ww. <- ww. <<- wSam # W[,i] # put in
#data1 <<- data.frame(data, ww, pweights)
allModels[[i]] <- gamlss(allFormula[[i]], weights=ww.*pweights., data = data, family = allFamily[[i]], control=g.control, ...)
W[,i] <- get.likelihood(allModels[[i]])
}
if(zero.component == TRUE) W[,KK] <- prob[KK]*ifelse(Y==0,1,0)
SumLik <- rowSums(W)
W <- W/SumLik
W <- ifelse(W<1e-10,0,W)
prob <- colSums(W)/N
# -2 * logLikelihood
newdv <- -2*sum(log(SumLik)) # the global deviance
olddv <- newdv + 1
iter.num <- 1
trace.print <- 0
dev.fits <- rep(0,control$n.cyc)
#-------------------------------------------------------------------------------
# if require model for prior propabilities pi's
modelPi <- if (pi.formula[[2]]==1) FALSE else TRUE
if (modelPi)
{
if (missing(data)) stop("the data argument is needed if pi's are modelled")
if(zero.component == TRUE)
{
dataKK <- expand.vc(data,KK)# expand data.
fac.fit <- gl(KK,N)
}
else
{
dataKK. <<- expand.vc(data,K)# expand data
fac.fit <- gl(K,N) # get the factor
}
res.var <- model.matrix(~fac.fit-1) # get the matrix of zeroes ans ones
form.prob <- update(pi.formula, res.var~.)# formula(paste("res.var~", pi.formula[[2]])) # get the formula
PROB <- matrix(prob, ncol=K, nrow=N, byrow=TRUE) # expand prior prob
}
## EM starts here---------------------------------------------------------------
while ( abs(olddv-newdv) > control$cc && iter.num < control$n.cyc ) # MS Wednesday, June 26, 2002
{
for (i in 1:K)
{ ww. <- ww. <<- W[,i]
allModels[[i]] <- gamlss(allFormula[[i]], weights=ww.*pweights., data = data, family = allFamily[[i]], control=g.control,...)
W[,i] <- if (modelPi) PROB[,i]*get.likelihood(allModels[[i]]) else prob[i]*get.likelihood(allModels[[i]])
}
if(zero.component == TRUE)
{
W[,KK] <- if (modelPi) PROB[,KK]*ifelse(Y==0,1,0) else prob[KK]*ifelse(Y==0,1,0)
}
SumLik <-rowSums(W)
W <- W/SumLik
W <- ifelse(W<1e-10,0,W)
di <- -2*log(SumLik)
olddv <- newdv
newdv <- sum(di)
if (modelPi)
{
dataKK.$wWw <- wWw <- as.vector(W)
dataKK.$res.var <- res.var
mult.mod <- multinom(form.prob, weights=wWw, data=dataKK., trace=FALSE)
fitted(mult.mod)[1:11,]
PROB <- fitted(mult.mod)[1:N,] # take only N of them
}
else { prob <- colSums(W)/N}
if (control$trace == TRUE) cat("GAMLSS-MX iteration ",iter.num, "Global deviance =", newdv, "\n" )
else if (is.numeric(control$trace))
{
trace.print <- trace.print+1
if (trace.print==control$trace)
{cat("GAMLSS-MX iteration ",iter.num, "Global deviance =", newdv, "\n" )
trace.print <- 0
}
}
dev.fits[iter.num] <- newdv
iter.num <- iter.num+1
}
if (control$plot==TRUE) plot(dev.fits[dev.fits!=0], type="l", xlab="EM iterations", ylab="-2 logLik")
## EM finish here---------------------------------------------------------------
## the output starts here
## residuals -------------------------------------------------------------------
WF <- WF2 <- matrix(0, ncol=K, nrow=N)
if (modelPi)
{
for (i in 1:K) WF[,i] <- PROB[,i]*get.the.p.function(allModels[[i]])
}
else
{
for (i in 1:K) WF[,i] <- prob[i]*get.the.p.function(allModels[[i]])
}
res <-qnorm(rowSums(WF))
for (i in 1:K) WF2[,i] <- W[,i]*get.the.p.function(allModels[[i]])
res2 <-qnorm(rowSums(WF2))
#-------------------------------------------------------------------------------
## degrees of freedom
df.fit <- 0
for (i in 1:K) df.fit <- df.fit + allModels[[i]]$df.fit
df.fit <- if(modelPi) df.fit+mult.mod$edf else df.fit+K-1 #
if (modelPi) colnames(PROB) <- paste("pi",seq(1,K), sep="")
familyAll <- rep("NULL",K)
for (i in 1:K) familyAll[i] <- allModels[[i]]$family[1]
## all output
out <- list()
out <- list(models = allModels,
model.pi = if (modelPi) mult.mod else NULL,
G.deviance = newdv,
df.fit = df.fit,
df.residual = N-df.fit,
post.prob = W,
prob = if (modelPi) PROB else prob,
family = familyAll,
call = gamlssMXcall,
aic = newdv+2*df.fit,
sbc = newdv+log(N)*df.fit,
K = K,
N = N,
weights = pweights.,
residuals = res,
Presiduals = res2,
seed = control$seed
)
class(out) <- list("gamlssMX", "gamlss")
if (modelPi) on.exit(rm(ww.,pweights.,dataKK., envir=sys.frame(0)))
else on.exit(rm(ww.,pweights., envir=sys.frame(0)))
out
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
MX.control <- function(cc=0.0001, n.cyc=200, trace=FALSE, seed=NULL, plot=TRUE, sample=NULL, ...)
{
##list(cc=0.0001, n.cyc=100, trace=TRUE),
## Control iteration for GAMLSS
## Mikis Stasinopoulos Monday, March 25, 2002 at 16:17
##
if(cc <= 0) {
warning("the value of cc supplied is zero or negative the default value of 0.0001 was used instead")
c.crit <- 0.0001}
if(n.cyc < 1) {
warning("the value of no cycles supplied is zero or negative the default value of 20 was used instead")
n.cyc <- 100}
if(is.logical(trace)) trace <- trace
else if (is.numeric(trace) & trace <= 0)
{warning("the value of trace supplied is less or equal t zero the default of 1 was used instead")
trace <- 1
}
seed <- if (is.null(seed) ) floor(runif(min=0, max=100000, n=1)) else seed
if (!is.null(sample))
{ if (sample>1||sample<0)
{
sample <-0.1
warning("the value of sample supplied is not probability the default of 0.10 was used instead")
}
else sample <- sample
}
plot <- if(is.logical(plot)) plot else TRUE
list(cc = cc, n.cyc = n.cyc, trace = trace, seed=seed, plot=plot, sample=sample)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# METHODS for "gamlssMX"
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
print.gamlssMX <- function (x, digits = max(3, getOption("digits") - 3), ...)
{
K <- x$K
cat("\nMixing Family: ", deparse(x$family), "\n", fill=TRUE)
cat("Fitting method: EM algorithm \n")
cat("\nCall: ", deparse(x$call,width.cutoff=50), "\n", fill=TRUE)
for (i in 1:K)
{
if ("mu" %in% x$models[[i]]$parameters)
{
cat("Mu Coefficients for model:", i, "\n")
print.default(format(coef(x$models[[i]], "mu"), digits = digits), print.gap = 2, quote = FALSE)
}
if ("sigma" %in% x$models[[i]]$parameters)
{
cat("Sigma Coefficients for model:", i, "\n")
print.default(format(coef(x$models[[i]], "sigma"), digits = digits), print.gap = 2, quote = FALSE)
}
if ("nu" %in% x$models[[i]]$parameters)
{
cat("Nu Coefficients for model:", i, "\n")
print.default(format(coef(x$models[[i]], "nu"), digits = digits), print.gap = 2, quote = FALSE)
}
if ("tau" %in% x$models[[i]]$parameters)
{
cat("Tau Coefficients for model:", i, "\n")
print.default(format(coef(x$models[[i]], "tau"), digits = digits), print.gap = 2, quote = FALSE)
}
}
if(is.matrix(x$prob))
{
cat("model for pi: \n")
print(coef(x$model.pi))
cat("\nEstimated probabilities: \n")
print(x$prob[1:3,])
cat("... \n" )
}
else cat("\nEstimated probabilities:", x$prob, "\n" )
cat("\nDegrees of Freedom for the fit:", x$df.fit, "Residual Deg. of Freedom ",
x$df.residual, "\n")
cat("Global Deviance: ", format(signif(x$G.deviance)),
"\n AIC: ", format(signif(x$aic)), "\n SBC: ",
format(signif(x$sbc)), "\n")
invisible(x)
}
################################################################################
# fitted.gamlssMX
################################################################################
# this gives the compoment fv if K is set otherwise average the componets
# using the fitted probabilities
# MS is this justified if mu not the mean?
# also does it makes sense for other parameters??
fitted.gamlssMX<-function (object, K=1, ... )
{
if (K%in%seq(1:object$K))
{
x <- fitted(object$models[[K]], ...)
}
else
{
WF <- matrix(0, ncol=object$K, nrow=object$N)
for (i in 1:object$K) WF[,i] <- object$prob[i]*fitted(object$models[[i]])
x <-rowSums(WF)
}
x
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# coefficients of the componets
coef.gamlssMX <- function(object, K=1, ...)
{
coef(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
formula.gamlssMX <- function(x, K=1, ...)
{
formula(x$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
model.matrix.gamlssMX <- function(object, K=1, ...)
{
model.matrix(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
terms.gamlssMX <- function(x, K=1, ...)
{
terms(x$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
predict.gamlssMX <- function(object, K=1, ...)
{
predict(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
residuals.gamlssMX <- function(object,type=c("prior","post"), ...)
{
type <- match.arg(type)
switch(type,"prior"=object$residuals,"post"= object$Presiduals)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this function fits n models with different starting values
# and select the one with smallest deviance
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this function fits n models with different starting values
# and select the one with smallest deviance
# last modification Thursday, April 30, 2009
gamlssMXfits <- function( n = 5,
formula = formula(data),
pi.formula = ~1,
family = "NO", # note is character
weights,
K = 2,
prob = NULL,
data = sys.parent(),
control = MX.control(),
g.control = gamlss.control(trace=FALSE),
zero.component = FALSE,
...)
{
gamlssMXfitscall <- match.call()
modellist <- list()
dev <- rep(0,length=n)
seed <- floor(runif(min=0,max=10000, n=n))
for (i in (1:n))
{
m1 <- try(gamlssMX(formula=formula, pi.formula = pi.formula, family = family ,K=K, prob=prob, data =
data, control=MX.control(seed=seed[i]), g.control = gamlss.control(trace=FALSE),
zero.component=zero.component)) #
if (any(class(m1)%in%"try-error"))
{
cat("model=", i, "failed", "\n")
dev[i] <- NA
modellist[[i]] <- NA
next
}
modellist[[i]] <- m1
dev[i] <- deviance(modellist[[i]])
cat("model=", i,"\n")
}
II <- which.min(dev)
model <- modellist[[II]]
gamlssMXfitscall$n <- NULL
gamlssMXfitscall[[1]] <- as.name("gamlssMX")
model$call <- gamlssMXfitscall
model$extra <- list(dev=dev, seed=seed, which=II)
model
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
update.gamlssMX <- function (object,
formula.,
...,
what = c("mu", "sigma", "nu", "tau"),
parameter= NULL,
evaluate = TRUE)
{
call <- object$call
if (is.null(call))
stop("need an object with call component")
extras <- match.call(expand.dots = FALSE)$...
if (!missing(formula.))
{
what <- if (!is.null(parameter)) {
match.arg(parameter, choices=c("mu", "sigma", "nu", "tau"))} else match.arg(what)
if (!what%in%object$par) stop(paste(what,"is not a parameter in the object","\n"))
if (what=="mu")
{ call$formula <- update.formula(formula(object,what), formula.) }
else
{
call[[paste(what,"formula",sep=".")]] <-
if (length(update.formula(formula(object,what), formula.))==2)
update.formula(formula(object,what), formula.)
else
update.formula(formula(object,what), formula.)[-2]
}
}
if (length(extras) > 0)
{
existing <- !is.na(match(names(extras), names(call)))
for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
if (any(!existing))
{
call <- c(as.list(call), extras[!existing])
call <- as.call(call)
}
}
if (evaluate)
eval(call, parent.frame())
else call
}
################################################################################
#
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this is a function to calculate the pdf of y fY(y)=p1*f1(y)+p2*f(y), of a MX model
# MS + BR
dMX <- function(y,
mu = list(mu1=1,mu2=5),
sigma = list(sigma1=1,sigma2=1),
nu = list(nu1=1,nu2=1),
tau = list(tau1=1,tau2=1),
pi = list(pi1=.2,pi2=.8),
# K = 2,
family = list(fam1="NO", fam2="NO"),
log = FALSE,
...)
{
# gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
# .gamlss.bi.list<-c("BI", "Binomial", "BB", "Beta Binomial")
## checking the probabilities
K <- length(pi)
sump <- sum(unlist(pi))
## set the length using the length of propabilities
# for (i in 1:K) sump <-pi[[i]]+ sump
if (!(9.9999999999>sump)&&(sump<1.00000000001)) stop(paste("the vector pi should sum to 1"))
## get the length of y and the length of the parameters
if(is.null(dim(y))) N <- length(y) else N <- dim(y)[1]
Prob <- matrix(1, nrow=N, ncol=K ) # rep(prob,rep(N,K))
for (i in 1:K) # go throught the components
{
fam <- as.gamlss.family(family[[i]])
fname <- fam$family[[1]]
DistPar <- fam$parameters
nopar <- length(DistPar)
dfun <- paste("d",fname,sep="")
if (fname%in%.gamlss.bi.list) # for binomial data
{
if (NCOL(y) == 1)
{
y <- y
bd <- rep(1, length(y))
}
else if (NCOL(y) == 2)
{
y <- y[,1]
bd <- y[,1]+y[,2]
}
else
{
stop("wrong response variable for binomial data")
}
}
switch(nopar,
{lik <- if (fname%in%.gamlss.bi.list){eval(call(dfun, x=y, bd=bd, mu=mu[[i]]))}
else eval(call(dfun, x=y, mu=mu[[i]]))
},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd, mu=mu[[i]], sigma=sigma[[i]]))
else eval(call(dfun, x=y, mu=mu[[i]], sigma=sigma[[i]]))
},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))
else eval(call(dfun,x=y, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))
else eval(call(dfun,x=y, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))})
Prob[,i]<-pi[[i]]*lik
} # finish go throught the components
fy <- rowSums(Prob)
fy <- if(log == FALSE) fy else log(fy)
fy
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this is a function to calculate the cdf y from a MX model
pMX <- function(q,
mu = list(mu1=1,mu2=5),
sigma = list(sigma1=1,sigma2=1),
nu = list(nu1=1,nu2=1),
tau = list(tau1=1,tau2=1),
pi = list(pi1=.2,pi2=.8),
# K = 2,
family = list(fam1="NO", fam2="NO"),
log = FALSE,
...)
{
# gamlss.bi.list<-c("BI", "Binomial", "BB", "Beta Binomial")
## checking the probabilities
K <- length(pi)
sump <- sum(unlist(pi))
if (sump!=1) stop(paste("the vector pi should sum to 1"))
# for (i in 1:K) sump <-pi[[i]]+ sump
if (!(9.9999999999>sump)&&(sump<1.00000000001)) stop(paste("the vector pi should sum to 1"))
## get the length of q and the length of the parameters
if(is.null(dim(q))) N <- length(q) else N <- dim(q)[1]
Prob <- matrix(1, nrow=N, ncol=K ) # rep(prob,rep(N,K))
for (i in 1:K)
{
fam <- as.gamlss.family(family[[i]])
fname <- fam$family[[1]]
DistPar <- fam$parameters
nopar <- length(DistPar)
dfun <- paste("p",fname,sep="")
if (fname%in%.gamlss.bi.list) # for binomial data
{
if (NCOL(q) == 1)
{
q <- q
bd <- rep(1, length(q))
}
else if (NCOL(q) == 2)
{
q <- q[,1]
bd <- q[,1]+q[,2]
}
else
{
stop("wrong response variable for binomial data")
}
}
switch(nopar,
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]], mu=mu[[i]]))
else eval(call(dfun, q=q, mu=mu[[i]]))
},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]], mu=mu[[i]], sigma=sigma[[i]]))
else eval(call(dfun, q=q, mu=mu[[i]], sigma=sigma[[i]]))
},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]], mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))
else eval(call(dfun,q=q, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]], mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))
else eval(call(dfun,q=q, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))})
Prob[,i]<-pi[[i]]*lik
}
fy <-rowSums(Prob)
fy <- if(log == FALSE) fy else log(fy)
fy
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# A function creating a function which can be use to plot a
# fitted mixture distrubution
# it needs a gamlssMX object and the observation number
getpdfMX <- function(object=NULL, observation=1)
{
if (class(object)[1]!="gamlssMX") stop("the object should be an gamlssMX object")
.gamlss.bi.list <- eval(quote(.gamlss.bi.list), envir = getNamespace("gamlss"))
K <- object$K
family <- object$family
ParamList <- list()
MU <- SIGMA <- NU<- TAU <- PI <- npar<- list()
modelForPi <- class(object$model.pi)[1]=="multinom"
for (i in 1:K)
{
ParamList[[i]] <- object$models[[i]]$parameters
if ("mu"%in%ParamList[[i]])
MU[[i]] <- fitted(object, K=i)[observation]
if ("sigma"%in%ParamList[[i]])
SIGMA[[i]] <- fitted(object, K=i, parameter="sigma")[observation]
if ("nu"%in%ParamList[[i]])
NU[[i]] <- fitted(object, K=i, parameter="nu")[observation]
if ("tau"%in%ParamList[[i]])
TAU[[i]] <- fitted(object, K=i, parameter="tau")[observation]
PI[[i]] <- if (modelForPi) object$prob[observation,i] else object$prob[i]
npar[[i]] <- length(ParamList[[i]])
}
npar <- unlist(npar)
if (any(npar==npar[1]))
nopar <- npar[1] else
stop("the different distributions should have same number of parameters for the function to work")
plotfun <- function(y)
{
if (object$family[1]%in%.gamlss.bi.list)
{
yy <- object$y
bd <- object$bd
switch(nopar,
dMX(y, bd=bd, mu=MU, pi=PI, family=family),
dMX(y, bd=bd, mu=MU, sigma=SIGMA, pi=PI, family=family),
dMX(y, bd=bd, mu=MU, sigma=SIGMA, nu=NU, pi=PI, family=family),
dMX(y, bd=bd, mu=MU, sigma=SIGMA, nu=NU, tau=TAU, pi=PI, family=family))
} else
{
switch(nopar,
dMX(y, mu=MU, pi=PI, family=family),
dMX(y, mu=MU, sigma=SIGMA, pi=PI, family=family),
dMX(y, mu=MU, sigma=SIGMA, nu=NU, pi=PI, family=family),
dMX(y, mu=MU, sigma=SIGMA, nu=NU, tau=TAU, pi=PI, family=family))
}
}
plotfun
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# at the moment only for -Inf +Inf continuous
# " discrete?
meanMX <- function(obj, observations=NULL, limit=c(-Inf, Inf))
{
# checking whether a gamlss object
if (class(obj)[1]!="gamlssMX") stop("the object should be an gamlssMX object")
#
obs <- if (is.null(observations)) 1:obj$N else observations
mean <- rep(0, length(obs))
j <- 1
for (i in obs)
{
fnLO <- getpdfMX(obj, observation=i)
f <-function(x) x*fnLO(x)
mean[j] <- integrate(f, limit[1], limit[2])$value
j <- j+1
}
mean
}
#
#
Try the gamlss.mx package in your browser
Any scripts or data that you put into this service are public.
gamlss.mx documentation built on May 29, 2024, 10:33 a.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 a general function for fitting mixures within gamlss
## created by Mikis Stasinopoulos and Bob Rigby
## LAST CHECKED 5-8-14
## latest change Tuesday, April 10, 2007 at 08:54
## Some of the models which can be fitted
## with the function gamlssNP can be identical
## In this version the prior probabilities pi can be modelled
## TO DO
## it needs a summary functions but the EM-algorithm SE's are not correct
## do we need EBP here?
## do we need means?
## the plotMX() need modification for BI type data Done but not checked
##
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
gamlssMX <- function ( formula = formula(data),
pi.formula = ~1,
family = "NO", # note is character
weights,
K = 2,
prob = NULL,
data = sys.parent(),
control = MX.control(...),
g.control = gamlss.control(trace=FALSE,...),
zero.component = FALSE,
...)
{
#-------------------------------------------------------------------------------
# this is to replicate rqres within gamlss enviroment DS Friday, March 31, 2006 at 10:30
rqres <- function (pfun = "pNO",
type = c("Continuous", "Discrete", "Mixed"),
censored = NULL,
ymin = NULL,
mass.p = NULL,
prob.mp = NULL,
y = y,
... )
{ }
body(rqres) <- eval(quote(body(rqres)), envir = getNamespace("gamlss"))
##------------------------------------------------------------------------------
##------------------------------------------------------------------------------
.gamlss.bi.list <- eval(quote(.gamlss.bi.list), envir = getNamespace("gamlss"))
# extra functions within -------------------------------------------------------
# for getting the commulative function------------------------------------------
get.the.p.function <- function(object, ...)
{
#gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
if (!is.gamlss(object)) stop(paste("This is not an gamlss object", "\n", ""))
fname <- object$family[[1]]
DistPar <- object$parameters
nopar <- length(DistPar)
dfun <- paste("p",fname,sep="")
# binomial denominators
switch(nopar,
{pfun <- if (fname%in%.gamlss.bi.list) eval(call(dfun,q=object$y, bd=object$bd, mu=fitted(object,"mu")))
else eval(call(dfun,q=object$y, mu=fitted(object,"mu")))},
{pfun <- if (fname%in%.gamlss.bi.list) eval(call(dfun,q=object$y, bd=object$bd, mu=fitted(object,"mu"),sigma=fitted(object,"sigma")))
else eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"))) },
{pfun <- if (fname%in%.gamlss.bi.list) eval(call(dfun,q=object$y, bd=object$bd, mu=fitted(object,"mu"),sigma=fitted(object,"sigma"), nu=fitted(object,"nu")))
else eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"), nu=fitted(object,"nu")))},
{pfun <- if (fname%in%.gamlss.bi.list) eval(call(dfun,q=object$y, bd=object$bd, mu=fitted(object,"mu"),sigma=fitted(object,"sigma"), nu=fitted(object,"nu"), tau=fitted(object,"tau")))
else eval(call(dfun,q=object$y, mu=fitted(object,"mu"), sigma=fitted(object,"sigma"), nu=fitted(object,"nu"), tau=fitted(object,"tau")))})
pfun
}
#-------------------------------------------------------------------------------
get.likelihood <- function(obj)
{
#gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
if (!is.gamlss(obj)) stop(paste("This is not an gamlss object", "\n", ""))
fname <- obj$family[[1]]
DistPar <- obj$parameters
nopar <- length(DistPar)
dfun <- paste("d",fname,sep="")
switch(nopar,
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun,x=obj$y, bd=obj$bd, mu=fitted(obj,"mu")))
else eval(call(dfun, x=obj$y, mu=fitted(obj)))},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun,x=obj$y, bd=obj$bd, mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma") ))
else eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"))) },
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun,x=obj$y, bd=obj$bd, mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu") ))
else eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu")))},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun,x=obj$y, bd=obj$bd, mu=fitted(obj,"mu"), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu"), tau=fitted(obj,"tau") ))
else eval(call(dfun,x=obj$y, mu=fitted(obj), sigma=fitted(obj,"sigma"), nu=fitted(obj,"nu"), tau=fitted(obj,"tau")))})
lik
}
#-------------------------------------------------------------------------------
expand.vc <- function(x,ni)
{
if (length(ni)==1)
{
if (ni<1) stop("ni should be greater than 1")
xx <- x
if (ni==1) return(xx)
for ( i in 2:ni) xx <- rbind(x,xx)
xx
}
else
{
n <- dim(x)[[1]]
c <- dim(x)[[2]]
xx <- NULL
for ( i in seq(1,n))
{
xx <- rbind(xx,matrix(rep(x[i,],ni[i]),ncol=c,byrow=TRUE))
}
xx
}
}
#-------------------------------------------------------------------------------
# the proper function starts here
#-------------------------------------------------------------------------------
#library(gamlss)
#library(nnet)
gamlssMXcall <- match.call() # the function call
## checking for NA in the data
if(!missing(data) & any(is.na(data)))
stop("The data contains NA's, use data = na.omit(mydata)")
## prob are the starting values of the pi's
if(zero.component == TRUE)
{ KK <- K+1
if (is.null(prob))
{
prob <- rep(1/KK, KK)
}
else
{
if (length(prob)!=KK) stop("the length of prob must be K+1(for zero component) ")
}
}
else
{
if (is.null(prob))
{
prob <- rep(1/K, K)
}
else
{
if (length(prob)!=K) stop("the length of prob must be K")
}
}
## expand the formula if not list
if (!is.list(formula))
{
allFormula <- vector("list",K)
for (i in 1:K) allFormula[[i]] <- formula
}
else allFormula <- formula
## expand the family if not list
if (!is.list(family))
{
allFamily <- vector("list",K)
for (i in 1:K) allFamily[[i]] <- family
}
else allFamily <- family
## check that the length for both the family and formula is K
if (length(allFamily)!=K) stop("the length of the list for the family is not equal the number of componets in the mixture")
if (length(allFormula)!=K) stop("the length of the list for the formula is not equal the number of componets in the mixture")
## get response for his length
Y <- model.extract(model.frame(allFormula[[1]],data=data),"response")
N <- if(is.null(dim(Y))) length(Y) else dim(Y)[1]# calculate the dimension for y
pweights. <- pweights. <<- if (missing(weights)) rep(1,N) else weights
## get things from control
set.seed(control$seed) # set the seed
prob.sample <- if (is.null(control$sample)) 10/N else control$sample
## creating the matrix of posterior probabilities (weights)
W <- matrix(1, nrow=N, ncol=K )
allModels <- vector("list",K)
## get starting values for the weights (posterior probabilities)
for (i in 1:K)
{ wSam <- sample(c(0,1), N, replace = TRUE, prob=c(1-min(.5,prob.sample),min(.5,prob.sample)))
ww. <- ww. <<- wSam # W[,i] # put in
#data1 <<- data.frame(data, ww, pweights)
allModels[[i]] <- gamlss(allFormula[[i]], weights=ww.*pweights., data = data, family = allFamily[[i]], control=g.control, ...)
W[,i] <- get.likelihood(allModels[[i]])
}
if(zero.component == TRUE) W[,KK] <- prob[KK]*ifelse(Y==0,1,0)
SumLik <- rowSums(W)
W <- W/SumLik
W <- ifelse(W<1e-10,0,W)
prob <- colSums(W)/N
# -2 * logLikelihood
newdv <- -2*sum(log(SumLik)) # the global deviance
olddv <- newdv + 1
iter.num <- 1
trace.print <- 0
dev.fits <- rep(0,control$n.cyc)
#-------------------------------------------------------------------------------
# if require model for prior propabilities pi's
modelPi <- if (pi.formula[[2]]==1) FALSE else TRUE
if (modelPi)
{
if (missing(data)) stop("the data argument is needed if pi's are modelled")
if(zero.component == TRUE)
{
dataKK <- expand.vc(data,KK)# expand data.
fac.fit <- gl(KK,N)
}
else
{
dataKK. <<- expand.vc(data,K)# expand data
fac.fit <- gl(K,N) # get the factor
}
res.var <- model.matrix(~fac.fit-1) # get the matrix of zeroes ans ones
form.prob <- update(pi.formula, res.var~.)# formula(paste("res.var~", pi.formula[[2]])) # get the formula
PROB <- matrix(prob, ncol=K, nrow=N, byrow=TRUE) # expand prior prob
}
## EM starts here---------------------------------------------------------------
while ( abs(olddv-newdv) > control$cc && iter.num < control$n.cyc ) # MS Wednesday, June 26, 2002
{
for (i in 1:K)
{ ww. <- ww. <<- W[,i]
allModels[[i]] <- gamlss(allFormula[[i]], weights=ww.*pweights., data = data, family = allFamily[[i]], control=g.control,...)
W[,i] <- if (modelPi) PROB[,i]*get.likelihood(allModels[[i]]) else prob[i]*get.likelihood(allModels[[i]])
}
if(zero.component == TRUE)
{
W[,KK] <- if (modelPi) PROB[,KK]*ifelse(Y==0,1,0) else prob[KK]*ifelse(Y==0,1,0)
}
SumLik <-rowSums(W)
W <- W/SumLik
W <- ifelse(W<1e-10,0,W)
di <- -2*log(SumLik)
olddv <- newdv
newdv <- sum(di)
if (modelPi)
{
dataKK.$wWw <- wWw <- as.vector(W)
dataKK.$res.var <- res.var
mult.mod <- multinom(form.prob, weights=wWw, data=dataKK., trace=FALSE)
fitted(mult.mod)[1:11,]
PROB <- fitted(mult.mod)[1:N,] # take only N of them
}
else { prob <- colSums(W)/N}
if (control$trace == TRUE) cat("GAMLSS-MX iteration ",iter.num, "Global deviance =", newdv, "\n" )
else if (is.numeric(control$trace))
{
trace.print <- trace.print+1
if (trace.print==control$trace)
{cat("GAMLSS-MX iteration ",iter.num, "Global deviance =", newdv, "\n" )
trace.print <- 0
}
}
dev.fits[iter.num] <- newdv
iter.num <- iter.num+1
}
if (control$plot==TRUE) plot(dev.fits[dev.fits!=0], type="l", xlab="EM iterations", ylab="-2 logLik")
## EM finish here---------------------------------------------------------------
## the output starts here
## residuals -------------------------------------------------------------------
WF <- WF2 <- matrix(0, ncol=K, nrow=N)
if (modelPi)
{
for (i in 1:K) WF[,i] <- PROB[,i]*get.the.p.function(allModels[[i]])
}
else
{
for (i in 1:K) WF[,i] <- prob[i]*get.the.p.function(allModels[[i]])
}
res <-qnorm(rowSums(WF))
for (i in 1:K) WF2[,i] <- W[,i]*get.the.p.function(allModels[[i]])
res2 <-qnorm(rowSums(WF2))
#-------------------------------------------------------------------------------
## degrees of freedom
df.fit <- 0
for (i in 1:K) df.fit <- df.fit + allModels[[i]]$df.fit
df.fit <- if(modelPi) df.fit+mult.mod$edf else df.fit+K-1 #
if (modelPi) colnames(PROB) <- paste("pi",seq(1,K), sep="")
familyAll <- rep("NULL",K)
for (i in 1:K) familyAll[i] <- allModels[[i]]$family[1]
## all output
out <- list()
out <- list(models = allModels,
model.pi = if (modelPi) mult.mod else NULL,
G.deviance = newdv,
df.fit = df.fit,
df.residual = N-df.fit,
post.prob = W,
prob = if (modelPi) PROB else prob,
family = familyAll,
call = gamlssMXcall,
aic = newdv+2*df.fit,
sbc = newdv+log(N)*df.fit,
K = K,
N = N,
weights = pweights.,
residuals = res,
Presiduals = res2,
seed = control$seed
)
class(out) <- list("gamlssMX", "gamlss")
if (modelPi) on.exit(rm(ww.,pweights.,dataKK., envir=sys.frame(0)))
else on.exit(rm(ww.,pweights., envir=sys.frame(0)))
out
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
MX.control <- function(cc=0.0001, n.cyc=200, trace=FALSE, seed=NULL, plot=TRUE, sample=NULL, ...)
{
##list(cc=0.0001, n.cyc=100, trace=TRUE),
## Control iteration for GAMLSS
## Mikis Stasinopoulos Monday, March 25, 2002 at 16:17
##
if(cc <= 0) {
warning("the value of cc supplied is zero or negative the default value of 0.0001 was used instead")
c.crit <- 0.0001}
if(n.cyc < 1) {
warning("the value of no cycles supplied is zero or negative the default value of 20 was used instead")
n.cyc <- 100}
if(is.logical(trace)) trace <- trace
else if (is.numeric(trace) & trace <= 0)
{warning("the value of trace supplied is less or equal t zero the default of 1 was used instead")
trace <- 1
}
seed <- if (is.null(seed) ) floor(runif(min=0, max=100000, n=1)) else seed
if (!is.null(sample))
{ if (sample>1||sample<0)
{
sample <-0.1
warning("the value of sample supplied is not probability the default of 0.10 was used instead")
}
else sample <- sample
}
plot <- if(is.logical(plot)) plot else TRUE
list(cc = cc, n.cyc = n.cyc, trace = trace, seed=seed, plot=plot, sample=sample)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# METHODS for "gamlssMX"
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
print.gamlssMX <- function (x, digits = max(3, getOption("digits") - 3), ...)
{
K <- x$K
cat("\nMixing Family: ", deparse(x$family), "\n", fill=TRUE)
cat("Fitting method: EM algorithm \n")
cat("\nCall: ", deparse(x$call,width.cutoff=50), "\n", fill=TRUE)
for (i in 1:K)
{
if ("mu" %in% x$models[[i]]$parameters)
{
cat("Mu Coefficients for model:", i, "\n")
print.default(format(coef(x$models[[i]], "mu"), digits = digits), print.gap = 2, quote = FALSE)
}
if ("sigma" %in% x$models[[i]]$parameters)
{
cat("Sigma Coefficients for model:", i, "\n")
print.default(format(coef(x$models[[i]], "sigma"), digits = digits), print.gap = 2, quote = FALSE)
}
if ("nu" %in% x$models[[i]]$parameters)
{
cat("Nu Coefficients for model:", i, "\n")
print.default(format(coef(x$models[[i]], "nu"), digits = digits), print.gap = 2, quote = FALSE)
}
if ("tau" %in% x$models[[i]]$parameters)
{
cat("Tau Coefficients for model:", i, "\n")
print.default(format(coef(x$models[[i]], "tau"), digits = digits), print.gap = 2, quote = FALSE)
}
}
if(is.matrix(x$prob))
{
cat("model for pi: \n")
print(coef(x$model.pi))
cat("\nEstimated probabilities: \n")
print(x$prob[1:3,])
cat("... \n" )
}
else cat("\nEstimated probabilities:", x$prob, "\n" )
cat("\nDegrees of Freedom for the fit:", x$df.fit, "Residual Deg. of Freedom ",
x$df.residual, "\n")
cat("Global Deviance: ", format(signif(x$G.deviance)),
"\n AIC: ", format(signif(x$aic)), "\n SBC: ",
format(signif(x$sbc)), "\n")
invisible(x)
}
################################################################################
# fitted.gamlssMX
################################################################################
# this gives the compoment fv if K is set otherwise average the componets
# using the fitted probabilities
# MS is this justified if mu not the mean?
# also does it makes sense for other parameters??
fitted.gamlssMX<-function (object, K=1, ... )
{
if (K%in%seq(1:object$K))
{
x <- fitted(object$models[[K]], ...)
}
else
{
WF <- matrix(0, ncol=object$K, nrow=object$N)
for (i in 1:object$K) WF[,i] <- object$prob[i]*fitted(object$models[[i]])
x <-rowSums(WF)
}
x
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# coefficients of the componets
coef.gamlssMX <- function(object, K=1, ...)
{
coef(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
formula.gamlssMX <- function(x, K=1, ...)
{
formula(x$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
model.matrix.gamlssMX <- function(object, K=1, ...)
{
model.matrix(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
terms.gamlssMX <- function(x, K=1, ...)
{
terms(x$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
predict.gamlssMX <- function(object, K=1, ...)
{
predict(object$models[[K]], ...)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
residuals.gamlssMX <- function(object,type=c("prior","post"), ...)
{
type <- match.arg(type)
switch(type,"prior"=object$residuals,"post"= object$Presiduals)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this function fits n models with different starting values
# and select the one with smallest deviance
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this function fits n models with different starting values
# and select the one with smallest deviance
# last modification Thursday, April 30, 2009
gamlssMXfits <- function( n = 5,
formula = formula(data),
pi.formula = ~1,
family = "NO", # note is character
weights,
K = 2,
prob = NULL,
data = sys.parent(),
control = MX.control(),
g.control = gamlss.control(trace=FALSE),
zero.component = FALSE,
...)
{
gamlssMXfitscall <- match.call()
modellist <- list()
dev <- rep(0,length=n)
seed <- floor(runif(min=0,max=10000, n=n))
for (i in (1:n))
{
m1 <- try(gamlssMX(formula=formula, pi.formula = pi.formula, family = family ,K=K, prob=prob, data =
data, control=MX.control(seed=seed[i]), g.control = gamlss.control(trace=FALSE),
zero.component=zero.component)) #
if (any(class(m1)%in%"try-error"))
{
cat("model=", i, "failed", "\n")
dev[i] <- NA
modellist[[i]] <- NA
next
}
modellist[[i]] <- m1
dev[i] <- deviance(modellist[[i]])
cat("model=", i,"\n")
}
II <- which.min(dev)
model <- modellist[[II]]
gamlssMXfitscall$n <- NULL
gamlssMXfitscall[[1]] <- as.name("gamlssMX")
model$call <- gamlssMXfitscall
model$extra <- list(dev=dev, seed=seed, which=II)
model
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
update.gamlssMX <- function (object,
formula.,
...,
what = c("mu", "sigma", "nu", "tau"),
parameter= NULL,
evaluate = TRUE)
{
call <- object$call
if (is.null(call))
stop("need an object with call component")
extras <- match.call(expand.dots = FALSE)$...
if (!missing(formula.))
{
what <- if (!is.null(parameter)) {
match.arg(parameter, choices=c("mu", "sigma", "nu", "tau"))} else match.arg(what)
if (!what%in%object$par) stop(paste(what,"is not a parameter in the object","\n"))
if (what=="mu")
{ call$formula <- update.formula(formula(object,what), formula.) }
else
{
call[[paste(what,"formula",sep=".")]] <-
if (length(update.formula(formula(object,what), formula.))==2)
update.formula(formula(object,what), formula.)
else
update.formula(formula(object,what), formula.)[-2]
}
}
if (length(extras) > 0)
{
existing <- !is.na(match(names(extras), names(call)))
for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
if (any(!existing))
{
call <- c(as.list(call), extras[!existing])
call <- as.call(call)
}
}
if (evaluate)
eval(call, parent.frame())
else call
}
################################################################################
#
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this is a function to calculate the pdf of y fY(y)=p1*f1(y)+p2*f(y), of a MX model
# MS + BR
dMX <- function(y,
mu = list(mu1=1,mu2=5),
sigma = list(sigma1=1,sigma2=1),
nu = list(nu1=1,nu2=1),
tau = list(tau1=1,tau2=1),
pi = list(pi1=.2,pi2=.8),
# K = 2,
family = list(fam1="NO", fam2="NO"),
log = FALSE,
...)
{
# gamlss.bi.list <- c("BI", "Binomial", "BB", "Beta Binomial")
# .gamlss.bi.list<-c("BI", "Binomial", "BB", "Beta Binomial")
## checking the probabilities
K <- length(pi)
sump <- sum(unlist(pi))
## set the length using the length of propabilities
# for (i in 1:K) sump <-pi[[i]]+ sump
if (!(9.9999999999>sump)&&(sump<1.00000000001)) stop(paste("the vector pi should sum to 1"))
## get the length of y and the length of the parameters
if(is.null(dim(y))) N <- length(y) else N <- dim(y)[1]
Prob <- matrix(1, nrow=N, ncol=K ) # rep(prob,rep(N,K))
for (i in 1:K) # go throught the components
{
fam <- as.gamlss.family(family[[i]])
fname <- fam$family[[1]]
DistPar <- fam$parameters
nopar <- length(DistPar)
dfun <- paste("d",fname,sep="")
if (fname%in%.gamlss.bi.list) # for binomial data
{
if (NCOL(y) == 1)
{
y <- y
bd <- rep(1, length(y))
}
else if (NCOL(y) == 2)
{
y <- y[,1]
bd <- y[,1]+y[,2]
}
else
{
stop("wrong response variable for binomial data")
}
}
switch(nopar,
{lik <- if (fname%in%.gamlss.bi.list){eval(call(dfun, x=y, bd=bd, mu=mu[[i]]))}
else eval(call(dfun, x=y, mu=mu[[i]]))
},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd, mu=mu[[i]], sigma=sigma[[i]]))
else eval(call(dfun, x=y, mu=mu[[i]], sigma=sigma[[i]]))
},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))
else eval(call(dfun,x=y, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, x=y, bd=bd, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))
else eval(call(dfun,x=y, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))})
Prob[,i]<-pi[[i]]*lik
} # finish go throught the components
fy <- rowSums(Prob)
fy <- if(log == FALSE) fy else log(fy)
fy
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# this is a function to calculate the cdf y from a MX model
pMX <- function(q,
mu = list(mu1=1,mu2=5),
sigma = list(sigma1=1,sigma2=1),
nu = list(nu1=1,nu2=1),
tau = list(tau1=1,tau2=1),
pi = list(pi1=.2,pi2=.8),
# K = 2,
family = list(fam1="NO", fam2="NO"),
log = FALSE,
...)
{
# gamlss.bi.list<-c("BI", "Binomial", "BB", "Beta Binomial")
## checking the probabilities
K <- length(pi)
sump <- sum(unlist(pi))
if (sump!=1) stop(paste("the vector pi should sum to 1"))
# for (i in 1:K) sump <-pi[[i]]+ sump
if (!(9.9999999999>sump)&&(sump<1.00000000001)) stop(paste("the vector pi should sum to 1"))
## get the length of q and the length of the parameters
if(is.null(dim(q))) N <- length(q) else N <- dim(q)[1]
Prob <- matrix(1, nrow=N, ncol=K ) # rep(prob,rep(N,K))
for (i in 1:K)
{
fam <- as.gamlss.family(family[[i]])
fname <- fam$family[[1]]
DistPar <- fam$parameters
nopar <- length(DistPar)
dfun <- paste("p",fname,sep="")
if (fname%in%.gamlss.bi.list) # for binomial data
{
if (NCOL(q) == 1)
{
q <- q
bd <- rep(1, length(q))
}
else if (NCOL(q) == 2)
{
q <- q[,1]
bd <- q[,1]+q[,2]
}
else
{
stop("wrong response variable for binomial data")
}
}
switch(nopar,
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]], mu=mu[[i]]))
else eval(call(dfun, q=q, mu=mu[[i]]))
},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]], mu=mu[[i]], sigma=sigma[[i]]))
else eval(call(dfun, q=q, mu=mu[[i]], sigma=sigma[[i]]))
},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]], mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))
else eval(call(dfun,q=q, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]]))},
{lik <- if (fname%in%.gamlss.bi.list) eval(call(dfun, q=q, bd=bd[[i]], mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))
else eval(call(dfun,q=q, mu=mu[[i]], sigma=sigma[[i]], nu=nu[[i]], tau=tau[[i]]))})
Prob[,i]<-pi[[i]]*lik
}
fy <-rowSums(Prob)
fy <- if(log == FALSE) fy else log(fy)
fy
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# A function creating a function which can be use to plot a
# fitted mixture distrubution
# it needs a gamlssMX object and the observation number
getpdfMX <- function(object=NULL, observation=1)
{
if (class(object)[1]!="gamlssMX") stop("the object should be an gamlssMX object")
.gamlss.bi.list <- eval(quote(.gamlss.bi.list), envir = getNamespace("gamlss"))
K <- object$K
family <- object$family
ParamList <- list()
MU <- SIGMA <- NU<- TAU <- PI <- npar<- list()
modelForPi <- class(object$model.pi)[1]=="multinom"
for (i in 1:K)
{
ParamList[[i]] <- object$models[[i]]$parameters
if ("mu"%in%ParamList[[i]])
MU[[i]] <- fitted(object, K=i)[observation]
if ("sigma"%in%ParamList[[i]])
SIGMA[[i]] <- fitted(object, K=i, parameter="sigma")[observation]
if ("nu"%in%ParamList[[i]])
NU[[i]] <- fitted(object, K=i, parameter="nu")[observation]
if ("tau"%in%ParamList[[i]])
TAU[[i]] <- fitted(object, K=i, parameter="tau")[observation]
PI[[i]] <- if (modelForPi) object$prob[observation,i] else object$prob[i]
npar[[i]] <- length(ParamList[[i]])
}
npar <- unlist(npar)
if (any(npar==npar[1]))
nopar <- npar[1] else
stop("the different distributions should have same number of parameters for the function to work")
plotfun <- function(y)
{
if (object$family[1]%in%.gamlss.bi.list)
{
yy <- object$y
bd <- object$bd
switch(nopar,
dMX(y, bd=bd, mu=MU, pi=PI, family=family),
dMX(y, bd=bd, mu=MU, sigma=SIGMA, pi=PI, family=family),
dMX(y, bd=bd, mu=MU, sigma=SIGMA, nu=NU, pi=PI, family=family),
dMX(y, bd=bd, mu=MU, sigma=SIGMA, nu=NU, tau=TAU, pi=PI, family=family))
} else
{
switch(nopar,
dMX(y, mu=MU, pi=PI, family=family),
dMX(y, mu=MU, sigma=SIGMA, pi=PI, family=family),
dMX(y, mu=MU, sigma=SIGMA, nu=NU, pi=PI, family=family),
dMX(y, mu=MU, sigma=SIGMA, nu=NU, tau=TAU, pi=PI, family=family))
}
}
plotfun
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# at the moment only for -Inf +Inf continuous
# " discrete?
meanMX <- function(obj, observations=NULL, limit=c(-Inf, Inf))
{
# checking whether a gamlss object
if (class(obj)[1]!="gamlssMX") stop("the object should be an gamlssMX object")
#
obs <- if (is.null(observations)) 1:obj$N else observations
mean <- rep(0, length(obs))
j <- 1
for (i in obs)
{
fnLO <- getpdfMX(obj, observation=i)
f <-function(x) x*fnLO(x)
mean[j] <- integrate(f, limit[1], limit[2])$value
j <- j+1
}
mean
}
#
#
Try the gamlss.mx 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.