Nothing
R/mle.R
In MAINT.Data: Model and Analyse Interval Data
Defines functions
ConfNextModel
ConfNestedBy
NextModel
NestedBy
MANOVAPermTest
getConfig
Documented in
ConfNestedBy
ConfNextModel
getConfig
MANOVAPermTest
NestedBy
NextModel
getConfig <- function(Config=NULL,...) Config
setMethod("mle",
signature(Sdt = "IData"),
# function(Sdt,Model=c("Normal","SKNormal","NrmandSKN"),CovCase=1:4,
# SelCrit=c("BIC","AIC"), k2max=1e6, OptCntrl=list(), ...)
function(Sdt, Model="Normal", CovCase="AllC", SelCrit=c("BIC","AIC"), k2max=1e6, OptCntrl=list(), ...)
{
# Model <- match.arg(Model)
if (Model!="Normal" && Model!="SKNormal" && Model!= "NrmandSKN")
stop("Wrong value for the 'Model' argument\n")
if (CovCase[1]=="AllC") CovCase <- 1:4
else for (Conf in CovCase) if (!is.element(Conf,1:4)) stop("Wrong value for the 'CovCase' argument\n")
SelCrit <- match.arg(SelCrit)
q <- Sdt@NIVar
p <- 2*q
n <- Sdt@NObs
limlnk2 <- log(k2max)
Config <- getConfig(...)
if (is.null(Config))
{
Config <- ifelse(CovCase==1,1,CovCase+1)
CovCaseArg <- TRUE
CovCaseMap <- c(1,NA,2,3,4)
} else {
CovCaseArg <- FALSE
CovCaseMap <- 1:5
}
if (Sdt@NIVar==1) {
CovCase <- q1CovCase(CovCase)
Config <- q1Config(Config)
}
if (Model!="SKNormal") {
Nres <- IdtNmle(Sdt,OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
}
if (Model!="Normal") {
# SNres <- IdtSNmle(Sdt,OptCntrl=OptCntrl,limlnk2=limlnk2,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
SNres <- IdtSNmle(Sdt,OptCntrl=OptCntrl,limlnk2=limlnk2,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,...)
}
if (Model=="Normal") { return(Nres) }
if (Model=="SKNormal") { return(SNres) }
if (Model=="NrmandSKN")
{
Xscld <- NULL
for (Conf in Config) {
CvCase <- CovCaseMap[Conf]
if (SNres@logLiks[CvCase] < Nres@logLiks[CvCase]) {
if (is.null(Xscld)) {
Xscld <- scale(cbind(Sdt@MidP,Sdt@LogR))
Xmean <- attr(Xscld,"scaled:center")
Xsd <- attr(Xscld,"scaled:scale")
XsdOutPrd <- outer(Xsd,Xsd)
lglikdif <- -n*sum(log(Xsd)) # difference between the log-likelihoods for the original and standardized data
}
# SNres@CovConfCases[[CvCase]] <- CorrSNSol(Nres,SNres,CvCase,Conf,Xscld,Xmean,Xsd,XsdOutPrd,lglikdif,limlnk2,OptCntrl)
SNres@CovConfCases[[CvCase]] <- CorrSNSol(Nres,SNres,CvCase,Conf,Xscld,Xmean,Xsd,XsdOutPrd,lglikdif,limlnk2,OptCntrl,...)
SNres@logLiks[CvCase] <- LogLik <- SNres@CovConfCases[[CvCase]]$logLik
nfreepar <- SKnpar(Conf,p,q)
SNres@AICs[CvCase] <- -2*LogLik + 2*nfreepar
SNres@BICs[CvCase] <- -2*LogLik + log(n)*nfreepar
}
}
AICs <- c(Nres@AICs,SNres@AICs)
BICs <- c(Nres@BICs,SNres@BICs)
if (SelCrit=="AIC") {
bestmod <- which.min(AICs)
} else {
if (SelCrit=="BIC") { bestmod <- which.min(BICs) }
}
return( new("IdtSngNandSNDE",NMod=Nres,SNMod=SNres,ModelNames=c(Nres@ModelNames,SNres@ModelNames),
ModelType=c(Nres@ModelType,SNres@ModelType),ModelConfig=c(Nres@ModelConfig,SNres@ModelConfig),
SelCrit=SelCrit,NIVar=Sdt@NIVar,logLiks=c(Nres@logLiks,SNres@logLiks),AICs=AICs,BICs=BICs,
BestModel=bestmod,SngD=TRUE)
)
}
}
)
setMethod("MANOVA",
signature(Sdt = "IData"),
function(Sdt, grouping, Model=c("Normal","SKNormal","NrmandSKN"), CovCase=1:4,
SelCrit=c("BIC","AIC"), Mxt=c("Hom","Het","Loc","Gen"), CVtol=1.0e-5, k2max=1e6,
OptCntrl=list(), onerror=c("stop","warning","silentNull"), ...)
{
onerror <- match.arg(onerror)
limlnk2 <- log(k2max)
if (!is.factor(grouping)) { return(error(onerror,"'grouping' is not a factor\n")) }
if ( Sdt@NObs != length(grouping)) {
return(error(onerror,"The numbers of data and partition observations are different\n"))
}
Mxt <- match.arg(Mxt)
Model <- match.arg(Model)
SelCrit <- match.arg(SelCrit)
H0res <- mle(Sdt,Model=Model,CovCase=CovCase,SelCrit=SelCrit,OptCntrl=OptCntrl,getvcov=FALSE,...)
Config <- getConfig(...)
if (is.null(Config))
{
Config <- ifelse(CovCase==1,1,CovCase+1)
CovCaseArg <- TRUE
CovCaseMap <- c(1,NA,2,3,4)
} else {
CovCaseArg <- FALSE
CovCaseMap <- 1:5
}
if (Sdt@NIVar==1) {
CovCase <- q1CovCase(CovCase)
Config <- q1Config(Config)
}
grouping <- factor(grouping,exclude=NULL)
nk <- as.numeric(table(grouping))
q <- Sdt@NIVar
p <- 2*q
n <- Sdt@NObs
k <- length(nk)
if (k==1) {
return( error(onerror,
"The data belongs to one single group. A partition into at least two different groups is required\n")
)
}
if (length(H0res@BestModel)==0) {
return( error(onerror,"Procedure MANOVA failed to find a valid null model mle estimate\n") )
}
if (Model=="Normal" && (Mxt=="Hom" || Mxt=="Loc") ) {
H1res <- IdtNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
} else if (Model=="SKNormal" && (Mxt=="Hom" || Mxt=="Loc") ) {
H1res <- IdtSNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
} else if (Model=="Normal" && (Mxt=="Het" || Mxt=="Gen") ) {
H1res <- IdtHetMxtNmle(Sdt,grouping,CVtol=CVtol,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
} else if (Model=="SKNormal" && (Mxt=="Het" || Mxt=="Gen") ) {
H1res <- IdtFDMxtSNmle(Sdt,grouping,CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
} else if (Model=="NrmandSKN") {
if (Mxt=="Hom" || Mxt=="Loc") {
SNH1res <- IdtSNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
NH1res <- IdtNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
SNH1res <- IdtSNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
} else if (Mxt=="Het" || Mxt=="Gen") {
NH1res <- IdtHetMxtNmle(Sdt,grouping,CVtol=CVtol,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
SNH1res <- IdtFDMxtSNmle(Sdt,grouping,CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
}
Xscld <- NULL
for (Conf in Config) {
CvCase <- CovCaseMap[Conf]
if (SNH1res@logLiks[CvCase] < NH1res@logLiks[CvCase]) {
if (is.null(Xscld)) {
Xscld <- scale(cbind(Sdt@MidP,Sdt@LogR))
Xmean <- attr(Xscld,"scaled:center")
Xsd <- attr(Xscld,"scaled:scale")
XsdOutPrd <- outer(Xsd,Xsd)
lglikdif <- -n*sum(log(Xsd)) # difference between the log-likelihoods for the original and standardized data
}
SNH1res@CovConfCases[[CvCase]] <-
CorrHetSNSol(NH1res,SNH1res,CvCase,Conf,Xscld,Xmean,Xsd,XsdOutPrd,grouping,Mxt,lglikdif,limlnk2,OptCntrl,getvcov=FALSE,...)
SNH1res@logLiks[CvCase] <- LogLik <- SNH1res@CovConfCases[[CvCase]]$logLik
nfreepar <- SKnpar(Conf,p,q)
SNH1res@AICs[CvCase] <- -2*LogLik + 2*nfreepar
SNH1res@BICs[CvCase] <- -2*LogLik + log(n)*nfreepar
}
}
AICs <- c(NH1res@AICs,SNH1res@AICs)
BICs <- c(NH1res@BICs,SNH1res@BICs)
if (SelCrit=="AIC") {
bestmod <- which.min(AICs)
} else if (SelCrit=="BIC") {
bestmod <- which.min(BICs)
}
H1res <- new("IdtMxNandSNDE",grouping=grouping,NMod=NH1res,SNMod=SNH1res,
ModelNames=c(NH1res@ModelNames,SNH1res@ModelNames),ModelType=c(NH1res@ModelType,SNH1res@ModelType),
ModelConfig=c(NH1res@ModelConfig,SNH1res@ModelConfig),SelCrit=SelCrit,NIVar=Sdt@NIVar,
logLiks=c(NH1res@logLiks,SNH1res@logLiks),AICs=AICs,BICs=BICs,BestModel=bestmod,SngD=FALSE,Ngrps=k)
}
if (length(H1res@BestModel)==0) {
return(error(onerror,"Procedure MANOVA failed to find a valid alternative model mle estimate\n"))
}
H0Ll <- H0res@logLiks[H1res@BestModel]
H1Ll <- H1res@logLiks[H1res@BestModel]
ChiSq <- 2*(H1Ll-H0Ll)
names(H0Ll) <- names(H1Ll) <- names(ChiSq) <- NULL
if (Mxt=="Hom" || Mxt=="Loc") {
df <- p*(k-1)
} else if (Mxt=="Het" || Mxt=="Gen") {
if ( Model=="Normal" || (Model=="NrmandSKN" && bestmod <= 5) )
{
BestConf <- H1res@ModelConfig[H1res@BestModel]
if (CovCaseArg && BestConf>1) BestConf <- BestConf+1 # Convert CovCase indices to Config indices
df <- npar(BestConf,p,q,Ngrps=k,Mxt="Het") - npar(BestConf,p,q,Ngrps=1,Mxt="Hom")
} else {
if (Model=="SKNormal") {
BestConf <- H1res@ModelConfig[H1res@BestModel]
} else {
BestConf <- SNH1res@ModelConfig[SNH1res@BestModel]
}
if (CovCaseArg && BestConf>1) BestConf <- BestConf+1 # Convert CovCase indices to Config indices
df <- SKnpar(BestConf,p,q,Ngrps=k,Mxt="GenMod") - SKnpar(BestConf,p,q,Ngrps=1)
}
}
pvalue <- pchisq(ChiSq,df,lower.tail=FALSE)
if (Model=="Normal" && (Mxt=="Hom" || Mxt=="Loc") ) {
return( new("IdtClMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="Normal" && (Mxt=="Het" || Mxt=="Gen") ) {
return( new("IdtHetNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="SKNormal" && (Mxt=="Hom" || Mxt=="Loc") ) {
return( new("IdtLocSNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="SKNormal" && (Mxt=="Het" || Mxt=="Gen") ) {
return( new("IdtGenSNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="NrmandSKN" && (Mxt=="Hom" || Mxt=="Loc") ) {
return( new("IdtLocNSNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="NrmandSKN" && (Mxt=="Het" || Mxt=="Gen") ) {
return( new("IdtGenNSNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
}
}
)
MANOVAPermTest <- function(MANOVAres, Sdt, grouping, nrep=200,
Model=c("Normal","SKNormal","NrmandSKN"), CovCase=1:4,
SelCrit=c("BIC","AIC"), Mxt=c("Hom","Het","Loc","Gen"), CVtol=1.0e-5, k2max=1e6,
OptCntrl=list(), onerror=c("stop","warning","silentNull"), ...)
{
if (!inherits(MANOVAres,"IdtMANOVA")) stop("Argument MANOVAres is not of class IdtMANOVA\n")
ChiSq <- MANOVAres@ChiSq
if (!inherits(Sdt,"IData")) stop("Argument Sdt is not of class IData\n")
if (!is.factor(grouping)) stop("Argument rouping is not a factor\n")
n <- Sdt@NObs
if (length(grouping) != n) stop("The numbers of data and partition observations are different\n")
Mxt <- match.arg(Mxt)
Model <- match.arg(Model)
SelCrit <- match.arg(SelCrit)
grouping <- factor(grouping,exclude=NULL)
grplvls <- levels(grouping)
nk <- as.numeric(table(grouping))
k <- length(nk)
cumsumnk <- cumsum(nk)
curgrouping <- grouping
empdist <- numeric(nrep)
for (rep in 1:nrep) {
permut <- sort.int(runif(n),index.return=TRUE)$ix
lstind <- 0
for (g in 1:k) {
frstind <- lstind+1
lstind <- cumsumnk[g]
curgrouping[permut[frstind:lstind]] <- grplvls[g]
}
empdist[rep] <- MANOVA(Sdt,curgrouping,Model,CovCase,SelCrit,Mxt,CVtol,k2max,OptCntrl,onerror,...)@ChiSq
}
pvalue <- length(which(ChiSq<empdist))/nrep
cat("Permutation p-value of MANOVA statistic",ChiSq,":",pvalue,"\n")
pvalue
}
setMethod("summary",
signature(object = "IdtE"),
function(object)
{
cat("Log likelihoods:\n")
print(object@logLiks)
if (object@SelCrit=="AIC") {
cat("Akaike Information Criteria:\n")
print(object@AICs)
} else if (object@SelCrit=="BIC") {
cat("Bayesian (Schwartz) Information Criteria:\n")
print(object@BICs)
}
cat("Selected model:\n")
print(names(object@BestModel))
invisible()
}
)
setMethod("show",
signature(object = "IdtE"),
function(object)
{
summary(object)
cat("\nSelected model parameter estimates:\n")
print(coef(object))
invisible()
}
)
setMethod("testMod",
signature(ModE = "IdtE"),
function(ModE, RestMod=ModE@ModelConfig[2]:length(ModE@ModelConfig), FullMod="Next")
{
if ( substr(ModE@ModelNames[1],2,8) == "ModCovC" || substr(ModE@ModelNames[1],3,9) == "ModCovC")
{
CovCaseArg <- TRUE
nCovCases <- 4
} else {
CovCaseArg <- FALSE
nCovCases <- 5
}
if (is.character(RestMod)) {
RestMod <- sapply(RestMod,function(Rmd) which(Rmd==ModE@ModelNames))
}
if (FullMod[1]!="Next" && FullMod[1]!="All" && is.character(FullMod)) {
FullMod <- sapply(FullMod,function(Fmd) which(Fmd==ModE@ModelNames))
}
if (is.element(ModE@ModelConfig[1],RestMod)) {
stop("Model",ModE@ModelNames,"can not be sepecified as a restricted model\n",
"since it is the most general model that has been estimated\n")
}
if (length(ModE@ModelType)>nCovCases) {
EType <- "NrmandSKN"
} else if (ModE@ModelType[1]=="Normal") {
EType <- "Normal"
} else if (ModE@ModelType[1]=="SkewNormal") {
EType <- "SKNormal"
}
if (is.character(RestMod)) {
RestMod <- sapply(RestMod,function(Rmd) which(Rmd==ModE@ModelNames))
}
if (FullMod[1]!="Next" && FullMod[1]!="All" && is.character(FullMod)) {
FullMod <- sapply(FullMod,function(Fmd) which(Fmd==ModE@ModelNames))
}
if (EType != "NrmandSKN" && is.element(ModE@ModelConfig[1],RestMod)) {
stop("Model",ModE@ModelNames[1],"can not be sepecified as a restricted model\n",
"since it is the most general model that has been estimated\n")
}
if (EType == "NrmandSKN" && is.element(nCovCases+1,RestMod)) {
MostGSNInd <- which(RestMod==nCovCases+1)
if (is.element(ModE@ModelConfig[1],RestMod)) {
RestMod <- RestMod[-MostGSNInd]
} else {
RestMod[MostGSNInd] <- ModE@ModelConfig[1]
}
}
TestRes <- list()
RestModels <- character()
FullModels <- character()
for (RMind in RestMod) {
if (FullMod[1]=="Next") {
FMindices <- NextModel(RMind,Model=EType,CovCaseArg=CovCaseArg)
} else if (FullMod[1]=="All") {
FMindices <- NestedBy(RMind,Model=EType,CovCaseArg=CovCaseArg)
} else {
FMindices <- intersect(FullMod,NestedBy(RMind,Model=EType,CovCaseArg=CovCaseArg))
}
for (FMind in FMindices) {
H0Ll <- ModE@logLiks[RMind]
H1Ll <- ModE@logLiks[FMind]
ChiSq <- 2*(H1Ll-H0Ll)
q <- ModE@NIVar
RConfig <- ModE@ModelConfig[RMind]
FConfig <- ModE@ModelConfig[FMind]
if (CovCaseArg) { # Convert CovCase indices to Config indices
if (RConfig>1) RConfig <- RConfig+1
if (FConfig>1) FConfig <- FConfig+1
}
if (ModE@SngD) {
if (ModE@ModelType[RMind]=="Normal") {
RMnpar <- npar(RConfig,2*q,q)
} else {
RMnpar <- SKnpar(RConfig,2*q,q)
}
if (ModE@ModelType[FMind]=="Normal") {
FMnpar <- npar(FConfig,2*q,q)
} else {
FMnpar <- SKnpar(FConfig,2*q,q)
}
} else {
if (ModE@ModelType[RMind]=="Normal") {
RMnpar <- npar(RConfig,2*q,q,Ngrps=ModE@Ngrps,Mxt="Het")
} else {
RMnpar <- SKnpar(RConfig,2*q,q,Ngrps=ModE@Ngrps,Mxt="GenMod")
}
if (ModE@ModelType[FMind]=="Normal") {
FMnpar <- npar(FConfig,2*q,q,Ngrps=ModE@Ngrps,Mxt="Het")
} else {
FMnpar <- SKnpar(FConfig,2*q,q,Ngrps=ModE@Ngrps,Mxt="GenMod")
}
}
df <- FMnpar - RMnpar
pvalue <- pchisq(ChiSq,df,lower.tail=FALSE)
resi <- new("LRTest",H0logLik=H0Ll,H1logLik=H1Ll,ChiSq=ChiSq,df=df,pvalue=pvalue)
TestRes <- c(TestRes,resi)
RestModels <- c(RestModels,ModE@ModelNames[RMind])
FullModels <- c(FullModels,ModE@ModelNames[FMind])
}
}
new("ConfTests",TestRes=TestRes,RestModels=RestModels,FullModels=FullModels)
}
)
NestedBy <- function(ModelInd,Model=c("Normal","SKNormal","NrmandSKN"),CovCaseArg)
{
Model <- match.arg(Model)
if (!CovCaseArg) { return(ConfNestedBy(ModelInd,Model)) }
if (Model=="Normal" || Model=="SKNormal")
{
if (ModelInd==1) { stop("Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==2 || ModelInd==3) { return(1) }
if (ModelInd==4) { return(1:3) }
} else if (Model=="NrmandSKN") {
if (ModelInd==1) { return(5) }
if (ModelInd==2) { return(c(1,5,6)) }
if (ModelInd==3) { return(c(1,5:7)) }
if (ModelInd==4) { return(c(1:3,5:8)) }
if (ModelInd==5) { stop("A Skew-Normal distribution with Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==6 || ModelInd==7) { return(5) }
if (ModelInd==8) { return(5:7) }
}
}
NextModel <- function(ModelInd,Model=c("Normal","SKNormal","NrmandSKN"),CovCaseArg)
{
Model <- match.arg(Model)
if (!CovCaseArg) { return(ConfNestedBy(ModelInd,Model)) }
if (Model=="Normal" || Model=="SKNormal")
{
if (ModelInd==1) { stop("Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==2 || ModelInd==3) { return(1) }
if (ModelInd==4) { return(2:3) }
} else if (Model=="NrmandSKN") {
if (ModelInd==1) { return(5) }
if (ModelInd==2) { return(c(1,6)) }
if (ModelInd==3) { return(c(1,7)) }
if (ModelInd==4) { return(c(2:3,8)) }
if (ModelInd==5) { stop("A Skew-Normal distribution with Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==6 || ModelInd==7) { return(5) }
if (ModelInd==8) { return(6:7) }
}
}
ConfNestedBy <- function(ModelInd,Model=c("Normal","SKNormal","NrmandSKN"))
{
Model <- match.arg(Model)
if (Model=="Normal" || Model=="SKNormal")
{
if (ModelInd==1) { stop("Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==2) { return(1) }
if (ModelInd==3 || ModelInd==4) { return(1:2) }
if (ModelInd==5) { return(1:4) }
} else if (Model=="NrmandSKN") {
if (ModelInd==1) { return(6) }
if (ModelInd==2) { return(c(1,6:7)) }
if (ModelInd==3) { return(c(1:2,6:8)) }
if (ModelInd==4) { return(c(1:2,6:7,9)) }
if (ModelInd==5) { return(c(1:4,6:10)) }
if (ModelInd==6) { stop("A Skew-Normal distribution with Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==7) { return(6) }
if (ModelInd==8 || ModelInd==9) { return(6:7) }
if (ModelInd==10) { return(6:9) }
}
}
ConfNextModel <- function(ModelInd,Model=c("Normal","SKNormal","NrmandSKN"))
{
Model <- match.arg(Model)
if (Model=="Normal" || Model=="SKNormal")
{
if (ModelInd==1) { stop("Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==2) { return(1) }
if (ModelInd==3 || ModelInd==4) { return(2) }
if (ModelInd==5) { return(3:4) }
} else if (Model=="NrmandSKN") {
if (ModelInd==1) { return(6) }
if (ModelInd==2) { return(c(1,7)) }
if (ModelInd==3) { return(c(2,8)) }
if (ModelInd==4) { return(c(2,9)) }
if (ModelInd==5) { return(c(3:4,10)) }
if (ModelInd==6) { stop("A Skew-Normal distribution with Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==7) { return(6) }
if (ModelInd==8 || ModelInd==9) { return(7) }
if (ModelInd==10) { return(8:9) }
}
}
setMethod("show",
signature(object = "ConfTests"),
function(object)
{
for (i in 1:length(object@TestRes)) {
if (!is.na(object@TestRes[[i]]@ChiSq))
{
cat("Testing Model",object@RestModels[i],"against alternative",object@FullModels[i],":\n")
print(object@TestRes[[i]])
}
}
invisible()
}
)
setMethod("BestModel",
signature(ModE = "IdtE"),
function(ModE,SelCrit=c("IdtCrt","BIC","AIC"))
{
SelCrit <- match.arg(SelCrit)
if (SelCrit == "IdtCrt") {
return(ModE@BestModel)
} else if (SelCrit == "AIC") {
return(which.min(ModE@AICs))
} else if (SelCrit == "BIC") {
return(which.min(ModE@BICs))
}
}
)
setMethod("coef",
signature(object = "IdtNDE"),
function(object,selmodel=BestModel(object),...)
{
# if (class(object)[1]=="IdtSngNDE" || class(object)[1]=="IdtMxNDE") {
if (inherits(object,"IdtSngNDE") || inherits(object,"IdtMxNDE")) {
return(list(mu=object@mleNmuE,Sigma=object@CovConfCases[[selmodel]]$mleSigE))
# } else if (class(object)[1]=="IdtSngNDRE" || class(object)[1]=="IdtMxNDRE") {
} else if (inherits(object,"IdtSngNDRE") || inherits(object,"IdtMxNDRE")) {
return(list(mu=object@RobNmuE,Sigma=object@CovConfCases[[selmodel]]$RobSigE))
}
}
)
setMethod("stdEr",
signature(x = "IdtNDE"),
function(x,selmodel=BestModel(x),...)
{
list(mu=x@mleNmuEse,Sigma=x@CovConfCases[[selmodel]]$mleSigEse)
}
)
setMethod("vcov",
signature(object = "IdtNDE"),
function(object,selmodel=BestModel(object),...)
{
object@CovConfCases[[selmodel]]$mlevcov
}
)
setMethod("mean", signature(x = "IdtNDE"), function(x) coef(x)$mu )
setMethod("var", signature(x ="IdtNDE"), function(x) coef(x)$Sigma )
setMethod("sd", signature(x ="IdtE"), function(x,na.rm=FALSE) sqrt(diag(var(x))) )
setMethod("AIC", signature(object="IdtE"),
function(object,...,k=2) {
if (missing(k)) nbarg <- 2
else nbarg <- 3
arguments <- match.call(expand.dots=TRUE)
if (length(arguments)==nbarg) {
aic <- object@AICs[object@BestModel]
if (k==2) return(aic)
else {
llik <- object@logLiks[object@BestModel]
return(-2*llik + (k/2)*(aic+2*llik))
}
}
nobjs <- length(arguments)-nbarg+1
aics <- numeric(nobjs)
objl <- vector("list",nobjs-1)
for (i in 1:nobjs) {
if (i==1) {
aics[1] <- object@AICs[object@BestModel]
names(aics)[1] <- names(object@AICs[object@BestModel])
} else {
objl[[i-1]] <- eval(arguments[[1+i]],sys.parent())
aics[i] <- objl[[i-1]]@AICs[objl[[i-1]]@BestModel]
names(aics)[i] <- names(objl[[i-1]]@AICs[objl[[i-1]]@BestModel])
}
}
if (k==2) return(aics)
else {
lliks <- numeric(nobjs)
for (i in 1:nobjs) {
if (i==1) lliks[1] <- object@logLiks[object@BestModel]
else lliks[i] <- objl[[i-1]]@logLiks[objl[[i-1]]@BestModel]
}
}
return(-2*lliks + (k/2)*(aics+2*lliks))
}
)
setMethod("BIC", signature(object="IdtE"),
function(object,...) {
arguments <- match.call(expand.dots=TRUE)
if (length(arguments)==2) return(object@BICs[object@BestModel])
nobjs <- length(arguments)-1
bics <- numeric(nobjs)
objl <- vector("list",nobjs-1)
for (i in 1:nobjs) {
if (i==1) {
bics[1] <- object@BICs[object@BestModel]
names(bics)[1] <- names(object@BICs[object@BestModel])
} else {
objl[[i-1]] <- eval(arguments[[1+i]],sys.parent())
bics[i] <- objl[[i-1]]@BICs[objl[[i-1]]@BestModel]
names(bics)[i] <- names(objl[[i-1]]@BICs[objl[[i-1]]@BestModel])
}
}
return(bics)
}
)
setMethod("logLik", signature(object="IdtE"),
function(object,...) {
arguments <- match.call(expand.dots=TRUE)
if (length(arguments)==2) return(object@logLiks[object@BestModel])
nobjs <- length(arguments)-1
logLiks <- numeric(nobjs)
objl <- vector("list",nobjs-1)
for (i in 1:nobjs) {
if (i==1) {
logLiks[1] <- object@logLiks[object@BestModel]
names(logLiks)[1] <- names(object@logLiks[object@BestModel])
} else {
objl[[i-1]] <- eval(arguments[[1+i]],sys.parent())
logLiks[i] <- objl[[i-1]]@logLiks[objl[[i-1]]@BestModel]
names(logLiks)[i] <- names(objl[[i-1]]@logLiks[objl[[i-1]]@BestModel])
}
}
return(logLiks)
}
)
setMethod("cor",
signature(x ="IdtNDE"),
function(x)
{
Sig <- coef(x)$Sigma
if (length(dim(Sig))==2) {
return(cov2cor(Sig))
} else if (length(dim(Sig))==3) {
return(array(apply(Sig,3,cov2cor),dim=dim(Sig),dimnames=dimnames(Sig)))
}
}
)
setMethod("vcov",
signature(object = "IdtMxNDE"),
function(object,selmodel=BestModel(object),group=NULL,...)
{
if (object@Hmcdt) {
object@CovConfCases[[selmodel]]$mlevcov
} else {
if (is.null(group))
{
warning(paste("vcov returned a three-dimensional array with a different var-cov matrix for each group,\n",
"which was identified by the level of the third array dimension\n"))
object@CovConfCases[[selmodel]]$mlevcov
} else {
object@CovConfCases[[selmodel]]$mlevcov[[,,group]]
}
}
}
)
setMethod("ObsLogLiks",
signature(object = "IdtSngNDE"),
function(object,Sdt,Conf=object@BestModel)
{
p <- 2*Sdt@NIVar
c0 <- -0.5*(p*log(2*pi))
Xdev <- scale(cbind(Sdt@MidP,Sdt@LogR),center=object@mleNmuE,scale=FALSE)
if (Conf!=5) {
SigISr <- t(backsolve(chol(object@CovConfCases[[Conf]]$mleSigE),diag(p)))
apply(Xdev,1,ILogLikNC1,SigmaSrInv=SigISr,const=c0+sum(log(diag(SigISr))))
} else {
IVar <- 1./diag(object@CovConfCases[[4]]$mleSigE)
apply(Xdev,1,ILogLikDNC,IVar=IVar,const=c0-0.5*sum(log(IVar)))
}
}
)
setMethod("show",
signature(object = "LRTest"),
function(object)
{
cat("Null Model log-likelihood:",object@H0logLik,"\n")
cat("Full Model log-likelihood:",object@H1logLik,"\n")
cat("Chi-squared statistic:",object@ChiSq,"\n")
cat("degrees of freedom:",object@df,"\n")
cat("p-value:",object@pvalue,"\n\n")
invisible()
}
)
setMethod("H1res", signature(object = "IdtMANOVA"), function(object) object@H1res)
setMethod("H0res", signature(object = "IdtMANOVA"), function(object) object@H0res)
setMethod("summary",
signature(object = "IdtMANOVA"),
function(object)
{
cat("Null Model Log likelihoods:\n")
print(object@H0res@logLiks)
cat("Full Model Log likelihoods:\n")
print(object@H1res@logLiks)
if (object@H1res@SelCrit=="AIC") {
cat("Full Model Akaike Information Criteria:\n")
print(object@H1res@AICs)
} else if (object@H1res@SelCrit=="BIC") {
cat("Full Model Bayesian (Schwartz) Information Criteria:\n")
print(object@H1res@BICs)
}
cat("Selected Model:\n")
print(names(object@H1res@BestModel))
cat("\n")
cat("Chi-squared statistic:",object@ChiSq,"\n")
cat("degrees of freedom:",object@df,"\n")
cat("p-value:",object@pvalue,"\n\n")
if ( length(object@grouping)<=30 ) {
cat("Note: Given the small sample size, the use of the Chi-square distribution may not be appropriate.\n",
"Alternatively, consider using the permutation test implemented in function MANOVAPermTest.\n",
"Note that this may take a long time.\n")
}
invisible()
}
)
setMethod("show",
signature(object = "IdtMANOVA"),
function(object) {
summary(object)
invisible()
}
)
Try the MAINT.Data package in your browser
Any scripts or data that you put into this service are public.
MAINT.Data documentation built on April 4, 2023, 9:09 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.
Defines functions ConfNextModel ConfNestedBy NextModel NestedBy MANOVAPermTest getConfig
Documented in ConfNestedBy ConfNextModel getConfig MANOVAPermTest NestedBy NextModel
getConfig <- function(Config=NULL,...) Config
setMethod("mle",
signature(Sdt = "IData"),
# function(Sdt,Model=c("Normal","SKNormal","NrmandSKN"),CovCase=1:4,
# SelCrit=c("BIC","AIC"), k2max=1e6, OptCntrl=list(), ...)
function(Sdt, Model="Normal", CovCase="AllC", SelCrit=c("BIC","AIC"), k2max=1e6, OptCntrl=list(), ...)
{
# Model <- match.arg(Model)
if (Model!="Normal" && Model!="SKNormal" && Model!= "NrmandSKN")
stop("Wrong value for the 'Model' argument\n")
if (CovCase[1]=="AllC") CovCase <- 1:4
else for (Conf in CovCase) if (!is.element(Conf,1:4)) stop("Wrong value for the 'CovCase' argument\n")
SelCrit <- match.arg(SelCrit)
q <- Sdt@NIVar
p <- 2*q
n <- Sdt@NObs
limlnk2 <- log(k2max)
Config <- getConfig(...)
if (is.null(Config))
{
Config <- ifelse(CovCase==1,1,CovCase+1)
CovCaseArg <- TRUE
CovCaseMap <- c(1,NA,2,3,4)
} else {
CovCaseArg <- FALSE
CovCaseMap <- 1:5
}
if (Sdt@NIVar==1) {
CovCase <- q1CovCase(CovCase)
Config <- q1Config(Config)
}
if (Model!="SKNormal") {
Nres <- IdtNmle(Sdt,OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
}
if (Model!="Normal") {
# SNres <- IdtSNmle(Sdt,OptCntrl=OptCntrl,limlnk2=limlnk2,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
SNres <- IdtSNmle(Sdt,OptCntrl=OptCntrl,limlnk2=limlnk2,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,...)
}
if (Model=="Normal") { return(Nres) }
if (Model=="SKNormal") { return(SNres) }
if (Model=="NrmandSKN")
{
Xscld <- NULL
for (Conf in Config) {
CvCase <- CovCaseMap[Conf]
if (SNres@logLiks[CvCase] < Nres@logLiks[CvCase]) {
if (is.null(Xscld)) {
Xscld <- scale(cbind(Sdt@MidP,Sdt@LogR))
Xmean <- attr(Xscld,"scaled:center")
Xsd <- attr(Xscld,"scaled:scale")
XsdOutPrd <- outer(Xsd,Xsd)
lglikdif <- -n*sum(log(Xsd)) # difference between the log-likelihoods for the original and standardized data
}
# SNres@CovConfCases[[CvCase]] <- CorrSNSol(Nres,SNres,CvCase,Conf,Xscld,Xmean,Xsd,XsdOutPrd,lglikdif,limlnk2,OptCntrl)
SNres@CovConfCases[[CvCase]] <- CorrSNSol(Nres,SNres,CvCase,Conf,Xscld,Xmean,Xsd,XsdOutPrd,lglikdif,limlnk2,OptCntrl,...)
SNres@logLiks[CvCase] <- LogLik <- SNres@CovConfCases[[CvCase]]$logLik
nfreepar <- SKnpar(Conf,p,q)
SNres@AICs[CvCase] <- -2*LogLik + 2*nfreepar
SNres@BICs[CvCase] <- -2*LogLik + log(n)*nfreepar
}
}
AICs <- c(Nres@AICs,SNres@AICs)
BICs <- c(Nres@BICs,SNres@BICs)
if (SelCrit=="AIC") {
bestmod <- which.min(AICs)
} else {
if (SelCrit=="BIC") { bestmod <- which.min(BICs) }
}
return( new("IdtSngNandSNDE",NMod=Nres,SNMod=SNres,ModelNames=c(Nres@ModelNames,SNres@ModelNames),
ModelType=c(Nres@ModelType,SNres@ModelType),ModelConfig=c(Nres@ModelConfig,SNres@ModelConfig),
SelCrit=SelCrit,NIVar=Sdt@NIVar,logLiks=c(Nres@logLiks,SNres@logLiks),AICs=AICs,BICs=BICs,
BestModel=bestmod,SngD=TRUE)
)
}
}
)
setMethod("MANOVA",
signature(Sdt = "IData"),
function(Sdt, grouping, Model=c("Normal","SKNormal","NrmandSKN"), CovCase=1:4,
SelCrit=c("BIC","AIC"), Mxt=c("Hom","Het","Loc","Gen"), CVtol=1.0e-5, k2max=1e6,
OptCntrl=list(), onerror=c("stop","warning","silentNull"), ...)
{
onerror <- match.arg(onerror)
limlnk2 <- log(k2max)
if (!is.factor(grouping)) { return(error(onerror,"'grouping' is not a factor\n")) }
if ( Sdt@NObs != length(grouping)) {
return(error(onerror,"The numbers of data and partition observations are different\n"))
}
Mxt <- match.arg(Mxt)
Model <- match.arg(Model)
SelCrit <- match.arg(SelCrit)
H0res <- mle(Sdt,Model=Model,CovCase=CovCase,SelCrit=SelCrit,OptCntrl=OptCntrl,getvcov=FALSE,...)
Config <- getConfig(...)
if (is.null(Config))
{
Config <- ifelse(CovCase==1,1,CovCase+1)
CovCaseArg <- TRUE
CovCaseMap <- c(1,NA,2,3,4)
} else {
CovCaseArg <- FALSE
CovCaseMap <- 1:5
}
if (Sdt@NIVar==1) {
CovCase <- q1CovCase(CovCase)
Config <- q1Config(Config)
}
grouping <- factor(grouping,exclude=NULL)
nk <- as.numeric(table(grouping))
q <- Sdt@NIVar
p <- 2*q
n <- Sdt@NObs
k <- length(nk)
if (k==1) {
return( error(onerror,
"The data belongs to one single group. A partition into at least two different groups is required\n")
)
}
if (length(H0res@BestModel)==0) {
return( error(onerror,"Procedure MANOVA failed to find a valid null model mle estimate\n") )
}
if (Model=="Normal" && (Mxt=="Hom" || Mxt=="Loc") ) {
H1res <- IdtNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
} else if (Model=="SKNormal" && (Mxt=="Hom" || Mxt=="Loc") ) {
H1res <- IdtSNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
} else if (Model=="Normal" && (Mxt=="Het" || Mxt=="Gen") ) {
H1res <- IdtHetMxtNmle(Sdt,grouping,CVtol=CVtol,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
} else if (Model=="SKNormal" && (Mxt=="Het" || Mxt=="Gen") ) {
H1res <- IdtFDMxtSNmle(Sdt,grouping,CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
} else if (Model=="NrmandSKN") {
if (Mxt=="Hom" || Mxt=="Loc") {
SNH1res <- IdtSNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
NH1res <- IdtNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
SNH1res <- IdtSNmle(Sdt,grouping,Type="HomMxt",CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
} else if (Mxt=="Het" || Mxt=="Gen") {
NH1res <- IdtHetMxtNmle(Sdt,grouping,CVtol=CVtol,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit)
SNH1res <- IdtFDMxtSNmle(Sdt,grouping,CVtol=CVtol,limlnk2=limlnk2,
OptCntrl=OptCntrl,CovCaseArg=CovCaseArg,Config=Config,SelCrit=SelCrit,getvcov=FALSE,...)
}
Xscld <- NULL
for (Conf in Config) {
CvCase <- CovCaseMap[Conf]
if (SNH1res@logLiks[CvCase] < NH1res@logLiks[CvCase]) {
if (is.null(Xscld)) {
Xscld <- scale(cbind(Sdt@MidP,Sdt@LogR))
Xmean <- attr(Xscld,"scaled:center")
Xsd <- attr(Xscld,"scaled:scale")
XsdOutPrd <- outer(Xsd,Xsd)
lglikdif <- -n*sum(log(Xsd)) # difference between the log-likelihoods for the original and standardized data
}
SNH1res@CovConfCases[[CvCase]] <-
CorrHetSNSol(NH1res,SNH1res,CvCase,Conf,Xscld,Xmean,Xsd,XsdOutPrd,grouping,Mxt,lglikdif,limlnk2,OptCntrl,getvcov=FALSE,...)
SNH1res@logLiks[CvCase] <- LogLik <- SNH1res@CovConfCases[[CvCase]]$logLik
nfreepar <- SKnpar(Conf,p,q)
SNH1res@AICs[CvCase] <- -2*LogLik + 2*nfreepar
SNH1res@BICs[CvCase] <- -2*LogLik + log(n)*nfreepar
}
}
AICs <- c(NH1res@AICs,SNH1res@AICs)
BICs <- c(NH1res@BICs,SNH1res@BICs)
if (SelCrit=="AIC") {
bestmod <- which.min(AICs)
} else if (SelCrit=="BIC") {
bestmod <- which.min(BICs)
}
H1res <- new("IdtMxNandSNDE",grouping=grouping,NMod=NH1res,SNMod=SNH1res,
ModelNames=c(NH1res@ModelNames,SNH1res@ModelNames),ModelType=c(NH1res@ModelType,SNH1res@ModelType),
ModelConfig=c(NH1res@ModelConfig,SNH1res@ModelConfig),SelCrit=SelCrit,NIVar=Sdt@NIVar,
logLiks=c(NH1res@logLiks,SNH1res@logLiks),AICs=AICs,BICs=BICs,BestModel=bestmod,SngD=FALSE,Ngrps=k)
}
if (length(H1res@BestModel)==0) {
return(error(onerror,"Procedure MANOVA failed to find a valid alternative model mle estimate\n"))
}
H0Ll <- H0res@logLiks[H1res@BestModel]
H1Ll <- H1res@logLiks[H1res@BestModel]
ChiSq <- 2*(H1Ll-H0Ll)
names(H0Ll) <- names(H1Ll) <- names(ChiSq) <- NULL
if (Mxt=="Hom" || Mxt=="Loc") {
df <- p*(k-1)
} else if (Mxt=="Het" || Mxt=="Gen") {
if ( Model=="Normal" || (Model=="NrmandSKN" && bestmod <= 5) )
{
BestConf <- H1res@ModelConfig[H1res@BestModel]
if (CovCaseArg && BestConf>1) BestConf <- BestConf+1 # Convert CovCase indices to Config indices
df <- npar(BestConf,p,q,Ngrps=k,Mxt="Het") - npar(BestConf,p,q,Ngrps=1,Mxt="Hom")
} else {
if (Model=="SKNormal") {
BestConf <- H1res@ModelConfig[H1res@BestModel]
} else {
BestConf <- SNH1res@ModelConfig[SNH1res@BestModel]
}
if (CovCaseArg && BestConf>1) BestConf <- BestConf+1 # Convert CovCase indices to Config indices
df <- SKnpar(BestConf,p,q,Ngrps=k,Mxt="GenMod") - SKnpar(BestConf,p,q,Ngrps=1)
}
}
pvalue <- pchisq(ChiSq,df,lower.tail=FALSE)
if (Model=="Normal" && (Mxt=="Hom" || Mxt=="Loc") ) {
return( new("IdtClMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="Normal" && (Mxt=="Het" || Mxt=="Gen") ) {
return( new("IdtHetNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="SKNormal" && (Mxt=="Hom" || Mxt=="Loc") ) {
return( new("IdtLocSNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="SKNormal" && (Mxt=="Het" || Mxt=="Gen") ) {
return( new("IdtGenSNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="NrmandSKN" && (Mxt=="Hom" || Mxt=="Loc") ) {
return( new("IdtLocNSNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
} else if (Model=="NrmandSKN" && (Mxt=="Het" || Mxt=="Gen") ) {
return( new("IdtGenNSNMANOVA",NIVar=Sdt@NIVar,grouping=grouping,H0res=H0res,H1res=H1res,
ChiSq=ChiSq,df=df,pvalue=pvalue,H0logLik=H0Ll,H1logLik=H1Ll) )
}
}
)
MANOVAPermTest <- function(MANOVAres, Sdt, grouping, nrep=200,
Model=c("Normal","SKNormal","NrmandSKN"), CovCase=1:4,
SelCrit=c("BIC","AIC"), Mxt=c("Hom","Het","Loc","Gen"), CVtol=1.0e-5, k2max=1e6,
OptCntrl=list(), onerror=c("stop","warning","silentNull"), ...)
{
if (!inherits(MANOVAres,"IdtMANOVA")) stop("Argument MANOVAres is not of class IdtMANOVA\n")
ChiSq <- MANOVAres@ChiSq
if (!inherits(Sdt,"IData")) stop("Argument Sdt is not of class IData\n")
if (!is.factor(grouping)) stop("Argument rouping is not a factor\n")
n <- Sdt@NObs
if (length(grouping) != n) stop("The numbers of data and partition observations are different\n")
Mxt <- match.arg(Mxt)
Model <- match.arg(Model)
SelCrit <- match.arg(SelCrit)
grouping <- factor(grouping,exclude=NULL)
grplvls <- levels(grouping)
nk <- as.numeric(table(grouping))
k <- length(nk)
cumsumnk <- cumsum(nk)
curgrouping <- grouping
empdist <- numeric(nrep)
for (rep in 1:nrep) {
permut <- sort.int(runif(n),index.return=TRUE)$ix
lstind <- 0
for (g in 1:k) {
frstind <- lstind+1
lstind <- cumsumnk[g]
curgrouping[permut[frstind:lstind]] <- grplvls[g]
}
empdist[rep] <- MANOVA(Sdt,curgrouping,Model,CovCase,SelCrit,Mxt,CVtol,k2max,OptCntrl,onerror,...)@ChiSq
}
pvalue <- length(which(ChiSq<empdist))/nrep
cat("Permutation p-value of MANOVA statistic",ChiSq,":",pvalue,"\n")
pvalue
}
setMethod("summary",
signature(object = "IdtE"),
function(object)
{
cat("Log likelihoods:\n")
print(object@logLiks)
if (object@SelCrit=="AIC") {
cat("Akaike Information Criteria:\n")
print(object@AICs)
} else if (object@SelCrit=="BIC") {
cat("Bayesian (Schwartz) Information Criteria:\n")
print(object@BICs)
}
cat("Selected model:\n")
print(names(object@BestModel))
invisible()
}
)
setMethod("show",
signature(object = "IdtE"),
function(object)
{
summary(object)
cat("\nSelected model parameter estimates:\n")
print(coef(object))
invisible()
}
)
setMethod("testMod",
signature(ModE = "IdtE"),
function(ModE, RestMod=ModE@ModelConfig[2]:length(ModE@ModelConfig), FullMod="Next")
{
if ( substr(ModE@ModelNames[1],2,8) == "ModCovC" || substr(ModE@ModelNames[1],3,9) == "ModCovC")
{
CovCaseArg <- TRUE
nCovCases <- 4
} else {
CovCaseArg <- FALSE
nCovCases <- 5
}
if (is.character(RestMod)) {
RestMod <- sapply(RestMod,function(Rmd) which(Rmd==ModE@ModelNames))
}
if (FullMod[1]!="Next" && FullMod[1]!="All" && is.character(FullMod)) {
FullMod <- sapply(FullMod,function(Fmd) which(Fmd==ModE@ModelNames))
}
if (is.element(ModE@ModelConfig[1],RestMod)) {
stop("Model",ModE@ModelNames,"can not be sepecified as a restricted model\n",
"since it is the most general model that has been estimated\n")
}
if (length(ModE@ModelType)>nCovCases) {
EType <- "NrmandSKN"
} else if (ModE@ModelType[1]=="Normal") {
EType <- "Normal"
} else if (ModE@ModelType[1]=="SkewNormal") {
EType <- "SKNormal"
}
if (is.character(RestMod)) {
RestMod <- sapply(RestMod,function(Rmd) which(Rmd==ModE@ModelNames))
}
if (FullMod[1]!="Next" && FullMod[1]!="All" && is.character(FullMod)) {
FullMod <- sapply(FullMod,function(Fmd) which(Fmd==ModE@ModelNames))
}
if (EType != "NrmandSKN" && is.element(ModE@ModelConfig[1],RestMod)) {
stop("Model",ModE@ModelNames[1],"can not be sepecified as a restricted model\n",
"since it is the most general model that has been estimated\n")
}
if (EType == "NrmandSKN" && is.element(nCovCases+1,RestMod)) {
MostGSNInd <- which(RestMod==nCovCases+1)
if (is.element(ModE@ModelConfig[1],RestMod)) {
RestMod <- RestMod[-MostGSNInd]
} else {
RestMod[MostGSNInd] <- ModE@ModelConfig[1]
}
}
TestRes <- list()
RestModels <- character()
FullModels <- character()
for (RMind in RestMod) {
if (FullMod[1]=="Next") {
FMindices <- NextModel(RMind,Model=EType,CovCaseArg=CovCaseArg)
} else if (FullMod[1]=="All") {
FMindices <- NestedBy(RMind,Model=EType,CovCaseArg=CovCaseArg)
} else {
FMindices <- intersect(FullMod,NestedBy(RMind,Model=EType,CovCaseArg=CovCaseArg))
}
for (FMind in FMindices) {
H0Ll <- ModE@logLiks[RMind]
H1Ll <- ModE@logLiks[FMind]
ChiSq <- 2*(H1Ll-H0Ll)
q <- ModE@NIVar
RConfig <- ModE@ModelConfig[RMind]
FConfig <- ModE@ModelConfig[FMind]
if (CovCaseArg) { # Convert CovCase indices to Config indices
if (RConfig>1) RConfig <- RConfig+1
if (FConfig>1) FConfig <- FConfig+1
}
if (ModE@SngD) {
if (ModE@ModelType[RMind]=="Normal") {
RMnpar <- npar(RConfig,2*q,q)
} else {
RMnpar <- SKnpar(RConfig,2*q,q)
}
if (ModE@ModelType[FMind]=="Normal") {
FMnpar <- npar(FConfig,2*q,q)
} else {
FMnpar <- SKnpar(FConfig,2*q,q)
}
} else {
if (ModE@ModelType[RMind]=="Normal") {
RMnpar <- npar(RConfig,2*q,q,Ngrps=ModE@Ngrps,Mxt="Het")
} else {
RMnpar <- SKnpar(RConfig,2*q,q,Ngrps=ModE@Ngrps,Mxt="GenMod")
}
if (ModE@ModelType[FMind]=="Normal") {
FMnpar <- npar(FConfig,2*q,q,Ngrps=ModE@Ngrps,Mxt="Het")
} else {
FMnpar <- SKnpar(FConfig,2*q,q,Ngrps=ModE@Ngrps,Mxt="GenMod")
}
}
df <- FMnpar - RMnpar
pvalue <- pchisq(ChiSq,df,lower.tail=FALSE)
resi <- new("LRTest",H0logLik=H0Ll,H1logLik=H1Ll,ChiSq=ChiSq,df=df,pvalue=pvalue)
TestRes <- c(TestRes,resi)
RestModels <- c(RestModels,ModE@ModelNames[RMind])
FullModels <- c(FullModels,ModE@ModelNames[FMind])
}
}
new("ConfTests",TestRes=TestRes,RestModels=RestModels,FullModels=FullModels)
}
)
NestedBy <- function(ModelInd,Model=c("Normal","SKNormal","NrmandSKN"),CovCaseArg)
{
Model <- match.arg(Model)
if (!CovCaseArg) { return(ConfNestedBy(ModelInd,Model)) }
if (Model=="Normal" || Model=="SKNormal")
{
if (ModelInd==1) { stop("Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==2 || ModelInd==3) { return(1) }
if (ModelInd==4) { return(1:3) }
} else if (Model=="NrmandSKN") {
if (ModelInd==1) { return(5) }
if (ModelInd==2) { return(c(1,5,6)) }
if (ModelInd==3) { return(c(1,5:7)) }
if (ModelInd==4) { return(c(1:3,5:8)) }
if (ModelInd==5) { stop("A Skew-Normal distribution with Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==6 || ModelInd==7) { return(5) }
if (ModelInd==8) { return(5:7) }
}
}
NextModel <- function(ModelInd,Model=c("Normal","SKNormal","NrmandSKN"),CovCaseArg)
{
Model <- match.arg(Model)
if (!CovCaseArg) { return(ConfNestedBy(ModelInd,Model)) }
if (Model=="Normal" || Model=="SKNormal")
{
if (ModelInd==1) { stop("Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==2 || ModelInd==3) { return(1) }
if (ModelInd==4) { return(2:3) }
} else if (Model=="NrmandSKN") {
if (ModelInd==1) { return(5) }
if (ModelInd==2) { return(c(1,6)) }
if (ModelInd==3) { return(c(1,7)) }
if (ModelInd==4) { return(c(2:3,8)) }
if (ModelInd==5) { stop("A Skew-Normal distribution with Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==6 || ModelInd==7) { return(5) }
if (ModelInd==8) { return(6:7) }
}
}
ConfNestedBy <- function(ModelInd,Model=c("Normal","SKNormal","NrmandSKN"))
{
Model <- match.arg(Model)
if (Model=="Normal" || Model=="SKNormal")
{
if (ModelInd==1) { stop("Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==2) { return(1) }
if (ModelInd==3 || ModelInd==4) { return(1:2) }
if (ModelInd==5) { return(1:4) }
} else if (Model=="NrmandSKN") {
if (ModelInd==1) { return(6) }
if (ModelInd==2) { return(c(1,6:7)) }
if (ModelInd==3) { return(c(1:2,6:8)) }
if (ModelInd==4) { return(c(1:2,6:7,9)) }
if (ModelInd==5) { return(c(1:4,6:10)) }
if (ModelInd==6) { stop("A Skew-Normal distribution with Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==7) { return(6) }
if (ModelInd==8 || ModelInd==9) { return(6:7) }
if (ModelInd==10) { return(6:9) }
}
}
ConfNextModel <- function(ModelInd,Model=c("Normal","SKNormal","NrmandSKN"))
{
Model <- match.arg(Model)
if (Model=="Normal" || Model=="SKNormal")
{
if (ModelInd==1) { stop("Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==2) { return(1) }
if (ModelInd==3 || ModelInd==4) { return(2) }
if (ModelInd==5) { return(3:4) }
} else if (Model=="NrmandSKN") {
if (ModelInd==1) { return(6) }
if (ModelInd==2) { return(c(1,7)) }
if (ModelInd==3) { return(c(2,8)) }
if (ModelInd==4) { return(c(2,9)) }
if (ModelInd==5) { return(c(3:4,10)) }
if (ModelInd==6) { stop("A Skew-Normal distribution with Configuration 1 is the most general model in this analysis\n") }
if (ModelInd==7) { return(6) }
if (ModelInd==8 || ModelInd==9) { return(7) }
if (ModelInd==10) { return(8:9) }
}
}
setMethod("show",
signature(object = "ConfTests"),
function(object)
{
for (i in 1:length(object@TestRes)) {
if (!is.na(object@TestRes[[i]]@ChiSq))
{
cat("Testing Model",object@RestModels[i],"against alternative",object@FullModels[i],":\n")
print(object@TestRes[[i]])
}
}
invisible()
}
)
setMethod("BestModel",
signature(ModE = "IdtE"),
function(ModE,SelCrit=c("IdtCrt","BIC","AIC"))
{
SelCrit <- match.arg(SelCrit)
if (SelCrit == "IdtCrt") {
return(ModE@BestModel)
} else if (SelCrit == "AIC") {
return(which.min(ModE@AICs))
} else if (SelCrit == "BIC") {
return(which.min(ModE@BICs))
}
}
)
setMethod("coef",
signature(object = "IdtNDE"),
function(object,selmodel=BestModel(object),...)
{
# if (class(object)[1]=="IdtSngNDE" || class(object)[1]=="IdtMxNDE") {
if (inherits(object,"IdtSngNDE") || inherits(object,"IdtMxNDE")) {
return(list(mu=object@mleNmuE,Sigma=object@CovConfCases[[selmodel]]$mleSigE))
# } else if (class(object)[1]=="IdtSngNDRE" || class(object)[1]=="IdtMxNDRE") {
} else if (inherits(object,"IdtSngNDRE") || inherits(object,"IdtMxNDRE")) {
return(list(mu=object@RobNmuE,Sigma=object@CovConfCases[[selmodel]]$RobSigE))
}
}
)
setMethod("stdEr",
signature(x = "IdtNDE"),
function(x,selmodel=BestModel(x),...)
{
list(mu=x@mleNmuEse,Sigma=x@CovConfCases[[selmodel]]$mleSigEse)
}
)
setMethod("vcov",
signature(object = "IdtNDE"),
function(object,selmodel=BestModel(object),...)
{
object@CovConfCases[[selmodel]]$mlevcov
}
)
setMethod("mean", signature(x = "IdtNDE"), function(x) coef(x)$mu )
setMethod("var", signature(x ="IdtNDE"), function(x) coef(x)$Sigma )
setMethod("sd", signature(x ="IdtE"), function(x,na.rm=FALSE) sqrt(diag(var(x))) )
setMethod("AIC", signature(object="IdtE"),
function(object,...,k=2) {
if (missing(k)) nbarg <- 2
else nbarg <- 3
arguments <- match.call(expand.dots=TRUE)
if (length(arguments)==nbarg) {
aic <- object@AICs[object@BestModel]
if (k==2) return(aic)
else {
llik <- object@logLiks[object@BestModel]
return(-2*llik + (k/2)*(aic+2*llik))
}
}
nobjs <- length(arguments)-nbarg+1
aics <- numeric(nobjs)
objl <- vector("list",nobjs-1)
for (i in 1:nobjs) {
if (i==1) {
aics[1] <- object@AICs[object@BestModel]
names(aics)[1] <- names(object@AICs[object@BestModel])
} else {
objl[[i-1]] <- eval(arguments[[1+i]],sys.parent())
aics[i] <- objl[[i-1]]@AICs[objl[[i-1]]@BestModel]
names(aics)[i] <- names(objl[[i-1]]@AICs[objl[[i-1]]@BestModel])
}
}
if (k==2) return(aics)
else {
lliks <- numeric(nobjs)
for (i in 1:nobjs) {
if (i==1) lliks[1] <- object@logLiks[object@BestModel]
else lliks[i] <- objl[[i-1]]@logLiks[objl[[i-1]]@BestModel]
}
}
return(-2*lliks + (k/2)*(aics+2*lliks))
}
)
setMethod("BIC", signature(object="IdtE"),
function(object,...) {
arguments <- match.call(expand.dots=TRUE)
if (length(arguments)==2) return(object@BICs[object@BestModel])
nobjs <- length(arguments)-1
bics <- numeric(nobjs)
objl <- vector("list",nobjs-1)
for (i in 1:nobjs) {
if (i==1) {
bics[1] <- object@BICs[object@BestModel]
names(bics)[1] <- names(object@BICs[object@BestModel])
} else {
objl[[i-1]] <- eval(arguments[[1+i]],sys.parent())
bics[i] <- objl[[i-1]]@BICs[objl[[i-1]]@BestModel]
names(bics)[i] <- names(objl[[i-1]]@BICs[objl[[i-1]]@BestModel])
}
}
return(bics)
}
)
setMethod("logLik", signature(object="IdtE"),
function(object,...) {
arguments <- match.call(expand.dots=TRUE)
if (length(arguments)==2) return(object@logLiks[object@BestModel])
nobjs <- length(arguments)-1
logLiks <- numeric(nobjs)
objl <- vector("list",nobjs-1)
for (i in 1:nobjs) {
if (i==1) {
logLiks[1] <- object@logLiks[object@BestModel]
names(logLiks)[1] <- names(object@logLiks[object@BestModel])
} else {
objl[[i-1]] <- eval(arguments[[1+i]],sys.parent())
logLiks[i] <- objl[[i-1]]@logLiks[objl[[i-1]]@BestModel]
names(logLiks)[i] <- names(objl[[i-1]]@logLiks[objl[[i-1]]@BestModel])
}
}
return(logLiks)
}
)
setMethod("cor",
signature(x ="IdtNDE"),
function(x)
{
Sig <- coef(x)$Sigma
if (length(dim(Sig))==2) {
return(cov2cor(Sig))
} else if (length(dim(Sig))==3) {
return(array(apply(Sig,3,cov2cor),dim=dim(Sig),dimnames=dimnames(Sig)))
}
}
)
setMethod("vcov",
signature(object = "IdtMxNDE"),
function(object,selmodel=BestModel(object),group=NULL,...)
{
if (object@Hmcdt) {
object@CovConfCases[[selmodel]]$mlevcov
} else {
if (is.null(group))
{
warning(paste("vcov returned a three-dimensional array with a different var-cov matrix for each group,\n",
"which was identified by the level of the third array dimension\n"))
object@CovConfCases[[selmodel]]$mlevcov
} else {
object@CovConfCases[[selmodel]]$mlevcov[[,,group]]
}
}
}
)
setMethod("ObsLogLiks",
signature(object = "IdtSngNDE"),
function(object,Sdt,Conf=object@BestModel)
{
p <- 2*Sdt@NIVar
c0 <- -0.5*(p*log(2*pi))
Xdev <- scale(cbind(Sdt@MidP,Sdt@LogR),center=object@mleNmuE,scale=FALSE)
if (Conf!=5) {
SigISr <- t(backsolve(chol(object@CovConfCases[[Conf]]$mleSigE),diag(p)))
apply(Xdev,1,ILogLikNC1,SigmaSrInv=SigISr,const=c0+sum(log(diag(SigISr))))
} else {
IVar <- 1./diag(object@CovConfCases[[4]]$mleSigE)
apply(Xdev,1,ILogLikDNC,IVar=IVar,const=c0-0.5*sum(log(IVar)))
}
}
)
setMethod("show",
signature(object = "LRTest"),
function(object)
{
cat("Null Model log-likelihood:",object@H0logLik,"\n")
cat("Full Model log-likelihood:",object@H1logLik,"\n")
cat("Chi-squared statistic:",object@ChiSq,"\n")
cat("degrees of freedom:",object@df,"\n")
cat("p-value:",object@pvalue,"\n\n")
invisible()
}
)
setMethod("H1res", signature(object = "IdtMANOVA"), function(object) object@H1res)
setMethod("H0res", signature(object = "IdtMANOVA"), function(object) object@H0res)
setMethod("summary",
signature(object = "IdtMANOVA"),
function(object)
{
cat("Null Model Log likelihoods:\n")
print(object@H0res@logLiks)
cat("Full Model Log likelihoods:\n")
print(object@H1res@logLiks)
if (object@H1res@SelCrit=="AIC") {
cat("Full Model Akaike Information Criteria:\n")
print(object@H1res@AICs)
} else if (object@H1res@SelCrit=="BIC") {
cat("Full Model Bayesian (Schwartz) Information Criteria:\n")
print(object@H1res@BICs)
}
cat("Selected Model:\n")
print(names(object@H1res@BestModel))
cat("\n")
cat("Chi-squared statistic:",object@ChiSq,"\n")
cat("degrees of freedom:",object@df,"\n")
cat("p-value:",object@pvalue,"\n\n")
if ( length(object@grouping)<=30 ) {
cat("Note: Given the small sample size, the use of the Chi-square distribution may not be appropriate.\n",
"Alternatively, consider using the permutation test implemented in function MANOVAPermTest.\n",
"Note that this may take a long time.\n")
}
invisible()
}
)
setMethod("show",
signature(object = "IdtMANOVA"),
function(object) {
summary(object)
invisible()
}
)
Try the MAINT.Data 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.