Nothing
context("MANOVA")
test_that("MANOVA works correctly for default homoscedastic Gaussian models", {
Agestrg <- substring(AbaloneIdt@ObsNames,first=3)
AbalClass <- factor(ifelse(Agestrg=="1-3"|Agestrg=="4-6"|Agestrg=="7-9","Young",
ifelse(Agestrg=="10-12"|Agestrg=="13-15"|Agestrg=="16-18","Adult","Old")
)
)
n <- nrow(AbaloneIdt)
q <- 2*ncol(AbaloneIdt) # Total number of MidPoints and LogRanges
k <- length(levels(AbalClass))
for (Cv in 1:4) {
AbMANOVA <- MANOVA(AbaloneIdt,AbalClass,CovCase=Cv)
AbMANOVAH0 <- H0res(AbMANOVA)
AbaIdtmle <- mle(AbaloneIdt,CovCase=Cv)
expect_equal(AbaIdtmle,AbMANOVAH0)
AbMANOVAH1 <- H1res(AbMANOVA)
SSCP <- matrix(0.,nrow=q,ncol=q)
rownames(SSCP) <- colnames(SSCP) <- c(names(MidPoints(AbaloneIdt)),names(LogRanges(AbaloneIdt)))
for (grp in levels(AbalClass)) {
grpIdt <- AbaloneIdt[AbalClass==grp,]
grpmle <- mle(grpIdt,CovCase=Cv)
expect_equal(mean(grpmle),mean(AbMANOVAH1)[grp,])
SSCP <- SSCP + nrow(grpIdt)*var(grpmle)
}
S <- SSCP/n
expect_equal(S,var(AbMANOVAH1))
}
} )
test_that("MANOVA computes correct standar errors for default homoscedastic Gaussian models", {
Agestrg <- substring(AbaloneIdt@ObsNames,first=3)
AbalClass <- factor(ifelse(Agestrg=="1-3"|Agestrg=="4-6"|Agestrg=="7-9","Young",
ifelse(Agestrg=="10-12"|Agestrg=="13-15"|Agestrg=="16-18","Adult","Old")
)
)
n <- nrow(AbaloneIdt)
q <- 2*ncol(AbaloneIdt) # Total number of MidPoints and LogRanges
k <- length(levels(AbalClass))
for (Cv in 1:4) {
AbMANOVA <- MANOVA(AbaloneIdt,AbalClass,CovCase=Cv)
AbMANOVAH0 <- H0res(AbMANOVA)
AbmeanStder <- sd(AbMANOVAH0) / sqrt(n)
expect_equal(stdEr(AbMANOVAH0)$mu,AbmeanStder)
vcovb_AbmeanStder <- sqrt(diag(vcov(AbMANOVAH0)[1:q,1:q]))
names(vcovb_AbmeanStder) <- names(AbmeanStder) <- NULL
expect_equal(vcovb_AbmeanStder,AbmeanStder)
mlecov <- var(AbMANOVAH0)
mlecov[mlecov==0.] <- NA
mlevar <- diag(mlecov)
AbcovStder <- sqrt( (mlecov^2 + outer(mlevar,mlevar)) / (n-1) )
expect_equal(stdEr(AbMANOVAH0)$Sigma,AbcovStder)
if (Cv==1) { # Implement later vcov checks for other configurations
vcovb_AbcovStder <- matrix(nrow=q,ncol=q)
vcovb_AbcovStder[lower.tri(vcovb_AbcovStder,diag=TRUE)] <- sqrt(diag(vcov(AbMANOVAH0)[-(1:q),-(1:q)]))
vcovb_AbcovStder[upper.tri(vcovb_AbcovStder)] <- t(vcovb_AbcovStder)[upper.tri(t(vcovb_AbcovStder))]
dimnames(vcovb_AbcovStder) <- dimnames(AbcovStder)
expect_equal(vcovb_AbcovStder,AbcovStder)
}
AbMANOVAH1 <- H1res(AbMANOVA)
nk <- as.numeric(table(AbalClass))
SSCP <- matrix(0.,nrow=q,ncol=q)
rownames(SSCP) <- colnames(SSCP) <- c(names(MidPoints(AbaloneIdt)),names(LogRanges(AbaloneIdt)))
for (grp in levels(AbalClass)) {
grpIdt <- AbaloneIdt[AbalClass==grp,]
grpmle <- mle(grpIdt,CovCase=Cv)
SSCP <- SSCP + nrow(grpIdt)*var(grpmle)
}
S <- SSCP/n
AbH1meannames <- list( levels(AbalClass), c(names(MidPoints(AbaloneIdt)),names(LogRanges(AbaloneIdt))) )
AbmeanStder <- matrix( rep(sqrt(diag(S)),each=k) / rep(sqrt(nk),q), nrow=k, ncol=q, dimnames=AbH1meannames )
expect_equal(AbmeanStder,stdEr(AbMANOVAH1)$mu)
vcovb_AbmeanStder <- matrix( sqrt(diag(vcov(AbMANOVAH1)[1:(k*q),1:(k*q)])), byrow=TRUE, nrow=k, ncol=q, dimnames=AbH1meannames )
expect_equal(vcovb_AbmeanStder,AbmeanStder)
mlecov <- var(AbMANOVAH1)
mlecov[mlecov==0.] <- NA
mlevar <- diag(mlecov)
AbcovStder <- sqrt( (mlecov^2 + outer(mlevar,mlevar)) / (n-k) )
expect_equal(stdEr(AbMANOVAH1)$Sigma,AbcovStder)
if (Cv==1) { # Implement later vcov checks for other configurations
vcovb_AbcovStder <- matrix(nrow=q,ncol=q)
vcovb_AbcovStder[lower.tri(vcovb_AbcovStder,diag=TRUE)] <- sqrt(diag(vcov(AbMANOVAH1)[-(1:(k*q)),-(1:(k*q))]))
vcovb_AbcovStder[upper.tri(vcovb_AbcovStder)] <- t(vcovb_AbcovStder)[upper.tri(t(vcovb_AbcovStder))]
dimnames(vcovb_AbcovStder) <- dimnames(AbcovStder)
expect_equal(vcovb_AbcovStder,AbcovStder)
}
}
} )
test_that("MANOVA works correctly for heteroscedastic Gaussian models", {
Agestrg <- substring(AbaloneIdt@ObsNames,first=3)
AbalClass <- factor(ifelse(Agestrg=="1-3"|Agestrg=="4-6"|Agestrg=="7-9","Young",
ifelse(Agestrg=="10-12"|Agestrg=="13-15"|Agestrg=="16-18","Adult","Old")
)
)
n <- nrow(AbaloneIdt)
q <- 2*ncol(AbaloneIdt) # Total number of MidPoints and LogRanges
k <- length(levels(AbalClass))
for (Cv in 1:4) {
AbMANOVA <- MANOVA(AbaloneIdt,AbalClass,Mxt="Het",CovCase=Cv)
AbMANOVAH0 <- H0res(AbMANOVA)
AbaIdtmle <- mle(AbaloneIdt,CovCase=Cv)
expect_equal(AbaIdtmle,AbMANOVAH0)
AbMANOVAH1 <- H1res(AbMANOVA)
for (grp in levels(AbalClass)) {
grpIdt <- AbaloneIdt[AbalClass==grp,]
grpmle <- mle(grpIdt,CovCase=Cv)
expect_equal(mean(grpmle),mean(AbMANOVAH1)[grp,])
expect_equal(var(grpmle),var(AbMANOVAH1)[,,grp])
}
}
} )
test_that("MANOVA computes correct standar errors for heteroscedastic Gaussian models", {
Agestrg <- substring(AbaloneIdt@ObsNames,first=3)
AbalClass <- factor(ifelse(Agestrg=="1-3"|Agestrg=="4-6"|Agestrg=="7-9","Young",
ifelse(Agestrg=="10-12"|Agestrg=="13-15"|Agestrg=="16-18","Adult","Old")
)
)
n <- nrow(AbaloneIdt)
q <- 2*ncol(AbaloneIdt) # Total number of MidPoints and LogRanges
k <- length(levels(AbalClass))
for (Cv in 1:4) {
AbMANOVA <- MANOVA(AbaloneIdt,AbalClass,Mxt="Het",CovCase=Cv)
AbMANOVAH0 <- H0res(AbMANOVA)
AbmeanStder <- sd(AbMANOVAH0) / sqrt(n)
expect_equal(stdEr(AbMANOVAH0)$mu,AbmeanStder)
vcovb_AbmeanStder <- sqrt(diag(vcov(AbMANOVAH0)[1:q,1:q]))
names(vcovb_AbmeanStder) <- names(AbmeanStder) <- NULL
expect_equal(vcovb_AbmeanStder,AbmeanStder)
mlecov <- var(AbMANOVAH0)
mlecov[mlecov==0.] <- NA
mlevar <- diag(mlecov)
AbcovStder <- sqrt( (mlecov^2 + outer(mlevar,mlevar)) / (n-1) )
expect_equal(stdEr(AbMANOVAH0)$Sigma,AbcovStder)
if (Cv==1) { # Implement later vcov checks for other configurations
vcovb_AbcovStder <- matrix(nrow=q,ncol=q)
vcovb_AbcovStder[lower.tri(vcovb_AbcovStder,diag=TRUE)] <- sqrt(diag(vcov(AbMANOVAH0)[-(1:q),-(1:q)]))
vcovb_AbcovStder[upper.tri(vcovb_AbcovStder)] <- t(vcovb_AbcovStder)[upper.tri(t(vcovb_AbcovStder))]
dimnames(vcovb_AbcovStder) <- dimnames(AbcovStder)
expect_equal(vcovb_AbcovStder,AbcovStder)
}
AbMANOVAH1 <- H1res(AbMANOVA)
nk <- as.numeric(table(AbalClass))
AbH1meannames <- list( levels(AbalClass), c(names(MidPoints(AbaloneIdt)),names(LogRanges(AbaloneIdt))) )
AbmeanStder <- matrix( nrow=k, ncol=q, dimnames=AbH1meannames )
for (g in 1:k) AbmeanStder[g,] <- sqrt( diag(var(AbMANOVAH1)[,,g]) / nk[g] )
expect_equal(AbmeanStder,stdEr(AbMANOVAH1)$mu)
vcovb_AbmeanStder <- matrix( nrow=k, ncol=q, dimnames=AbH1meannames )
for (g in 1:k) vcovb_AbmeanStder[g,] <- sqrt( diag(suppressWarnings(vcov(AbMANOVAH1))[1:q,1:q,g]) )
expect_equal(vcovb_AbmeanStder,AbmeanStder)
mlecov <- var(AbMANOVAH1)
mlecov[mlecov==0.] <- NA
for (g in 1:k) {
mlevar <- diag(mlecov[,,g])
AbcovStder <- sqrt( (mlecov[,,g]^2 + outer(mlevar,mlevar)) / (nk[g]-1) )
expect_equal(stdEr(AbMANOVAH1)$Sigma[,,g],AbcovStder)
if (Cv==1) { # Implement later vcov checks for other configurations
vcovb_AbcovStder <- matrix(nrow=q,ncol=q)
vcovb_AbcovStder[lower.tri(vcovb_AbcovStder,diag=TRUE)] <- sqrt( diag(suppressWarnings(vcov(AbMANOVAH1))[-(1:q),-(1:q),g]) )
vcovb_AbcovStder[upper.tri(vcovb_AbcovStder)] <- t(vcovb_AbcovStder)[upper.tri(t(vcovb_AbcovStder))]
dimnames(vcovb_AbcovStder) <- dimnames(AbcovStder)
expect_equal(vcovb_AbcovStder,AbcovStder)
}
}
}
} )
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.