Nothing
##############################################################################################
##################### F-MANOVA USING THE SGD APPROXIMATION FOR THE PRODUCT ###################
##############################################################################################
#' Computes a Mult-FANOVA model by an approximation
#' @param formula a description of the model to be fitted.
#' @param dataset the data frame containing all the variables of the model.
#' @param data.fuzzified the fuzzified data set constructed by a call to the function FUZZ or the function GFUZZ, or a similar matrix.
#' @param sig a numerical value representing the significance level of the test.
#' @param breakpoints a positive arbitrary integer representing the number of breaks chosen to build the numerical alpha-cuts. It is fixed to 100 by default.
#' @param index.var the column index of the considered variable for which the output will be printed. It is an argument of the Mult-FANOVA models by the exact and the approximation methods only.
#' @param int.method the method of numerical integration. It is set by default to the Simpson method, i.e. int.method="int.simpson".
#' @param plot fixed by default to "TRUE". plot="FALSE" if a plot of the fuzzy number is not required.
#' @return Returns a list of all the arguments of the function, the total, treatment and residuals sums of squares, the coefficients of the model, the test statistics with the corresponding p-values, and the decision made.
#' @importFrom FuzzyNumbers core
#' @importFrom stats model.frame
#' @importFrom stats model.matrix
#' @importFrom stats fitted.values
#' @importFrom stats model.response
#' @importFrom stats complete.cases
#' @importFrom stats coefficients
#' @importFrom FuzzyNumbers PiecewiseLinearFuzzyNumber
#' @importFrom graphics lines
#' @importFrom graphics legend
# #' @export
FMANOVA.approximation <- function(formula, dataset, data.fuzzified, sig=0.05, breakpoints= 100, index.var =NA, int.method = "int.simpson", plot = TRUE){
# START OF THE ALGORITHM
########################
if(is.trfuzzification(data.fuzzified) == TRUE){data.fuzzified <- tr.gfuzz(data.fuzzified, breakpoints = breakpoints)}
if(is.fuzzification(data.fuzzified) == FALSE){stop("Problems with the fuzzification matrix")}
breakpoints <- ncol(data.fuzzified) - 1
v <- c("TrapezoidalFuzzyNumber", "PowerFuzzyNumber", "PiecewiseLinearFuzzyNumber", "DiscontinuousFuzzyNumber", "FuzzyNumber")
if (unique(class(sig) %in% v) == TRUE){sig <- core(sig)[1]} else if (is.alphacuts(sig) == TRUE){sig <- sig[nrow(sig),1]
} else if (is.na(sig) == TRUE){stop("Significance level not defined")}
# Initialization of the datasets
mf <- model.frame(formula, dataset)
if (ncol(mf) == 2){stop("The FANOVA function should be used for this model")}
if (length(which((lapply(mf, nlevels)[1:ncol(mf)] > 2) == TRUE)) == 0){
data <- as.data.frame(model.matrix(mf, dataset))
} else {
dataset[,] <- lapply(dataset[,], as.numeric)
mf <- model.frame(formula, dataset)
data <- as.data.frame(model.matrix(mf, dataset))
}
Yc <- as.matrix(model.response(mf))
ok <- complete.cases(data, Yc)
data <- data[ok,]
Y <- data.fuzzified
data[,] <- lapply(data[,], factor)
if (colnames(data)[1] != "(Intercept)"){nc = ncol(data)} else if (colnames(data)[1] == "(Intercept)") {nc = ncol(data) - 1}
r <- matrix(rep(0), ncol=1, nrow = nc)
for(u in 2:ncol(data)){r[u-1,1] <- nlevels(data[,u])
#data[,u] <- relevel(data[,u], ref = r[u-1,1]) #2
# lapply(data[,], nlevels)[-1]
}
nt<- matrix(rep(0), ncol=1, nrow=nc)
nt[,] <- nrow(data[,])
ni <- matrix(rep(0), nrow= nc, ncol = max(r))
for(u in 2:(ncol(data))){ni[u-1,] <- table(data[,u])} #1:r[u-1,1]] <- table(data[,u])}
b <- breakpoints+1
mat.means <- array(rep(0), dim=c(nc,max(r)*b,2))
dist.mat.means <- matrix(rep(0), nrow = nc, ncol= max(r))
Y.. <- array(rep(0), dim=c(nc,b,2))
for (u in 2:ncol(data)){
WMS <- 0
for(v in 1:r[u-1,1]){
Part.mean <- Fuzzy.sample.mean(Y[which(data[,u]==levels(data[,u])[v]),,])
mat.means[u-1,((b*(v-1)+1):(b*(v-1)+b)),1] <- Part.mean[,1]
mat.means[u-1,((b*(v-1)+1):(b*(v-1)+b)),2] <- rev(Part.mean[,2])
#dist.mat.means[u-1,v] <- distance(Part.mean, TrapezoidalFuzzyNumber(0,0,0,0), type = distance.type, i=i, j=j, theta = theta, p=p, q=q, breakpoints = breakpoints)
LY.. <- Part.mean*ni[u-1,v]
WMS <- WMS + LY../nt[u-1]
}
Y..[u-1,,1] <- WMS[,1]
Y..[u-1,,2] <- rev(WMS[,2])
}
df <- matrix(rep(0), nrow=nc, ncol =1)
#E <- matrix(rep(0), nrow=nc, ncol =4)
#T <- matrix(rep(0), nrow=nc, ncol =4)
#H <- matrix(rep(0), nrow=nc, ncol =4)
E <- array(rep(0), dim=c(nc,b,2))
#T <- matrix(rep(0), nrow = nc, ncol=1)
T <- array(rep(0), dim=c(nc,b,2))
#H<- matrix(rep(0), nrow = nc, ncol=1)
H <- array(rep(0), dim=c(nc,b,2))
for(u in 2:(ncol(data))){
SE <- 0
ST <- 0
SH <- 0
E11S <- 0
T11S <- 0
H11S <- 0
T11.mean <- cbind(Y..[u-1,,1], rev(Y..[u-1,,2]))
colnames(T11.mean) <- c("L","U")
for(v in 1:r[u-1,1]){
Group <- Y[which(data[,u] == levels(data[,u])[v]),,]
E11 <- array(rep(0), dim=c(nrow(Group),b,2))
T11 <- array(rep(0), dim=c(nrow(Group),b,2))
H11 <- array(rep(0), dim=c(nc,b,2))
E11.mean <- cbind(mat.means[u-1,((b*(v-1)+1)):((b*(v-1)+b)),1], rev(mat.means[u-1,((b*(v-1)+1)):((b*(v-1)+b)),2]))
colnames(E11.mean) <- c("L","U")
for (w in 1:nrow(Group)){
E11.part <- cbind(Group[w,,1], rev(Group[w,,2]))
colnames(E11.part) <- c("L","U")
E11 <- Fuzzy.Difference(E11.part, E11.mean, alphacuts=TRUE, breakpoints = breakpoints)
TE11 <- c(E11[1,1], E11[(breakpoints+1),1], E11[(breakpoints+1),2], E11[1,2])
if (is.unsorted(TE11) == TRUE){TE11 <- sort(TE11)}
E11S <- E11S + (PiecewiseLinearFuzzyNumber(TE11[1], TE11[2], TE11[3], TE11[4]))*(PiecewiseLinearFuzzyNumber(TE11[1], TE11[2], TE11[3], TE11[4]))
T11 <- Fuzzy.Difference(E11.part, T11.mean, alphacuts=TRUE, breakpoints = breakpoints)
TT11 <- c(T11[1,1], T11[(breakpoints+1),1], T11[(breakpoints+1),2], T11[1,2])
if (is.unsorted(TT11) == TRUE){TT11 <- sort(TT11)}
T11S <- T11S + (PiecewiseLinearFuzzyNumber(TT11[1], TT11[2], TT11[3], TT11[4]))*(PiecewiseLinearFuzzyNumber(TT11[1], TT11[2], TT11[3], TT11[4]))
}
H11 <- Fuzzy.Difference(E11.mean, T11.mean, alphacuts=TRUE, breakpoints = breakpoints)
TH11 <- c(H11[1,1], H11[(breakpoints+1),1], H11[(breakpoints+1),2], H11[1,2])
if (is.unsorted(TH11) == TRUE){TH11 <- sort(TH11)}
H11S <- ((PiecewiseLinearFuzzyNumber(TH11[1], TH11[2], TH11[3], TH11[4]))*(PiecewiseLinearFuzzyNumber(TH11[1], TH11[2], TH11[3], TH11[4])))*ni[u-1,v]
SE = SE + E11S
ST = ST + T11S
SH = SH + H11S
H11 <- NULL
T11 <- NULL
E11 <- NULL
}
E[u-1,,1] <- alphacut(SE, seq(0,1,1/breakpoints))[,1]
E[u-1,,2] <- rev(alphacut(SE, seq(0,1,1/breakpoints))[,2])
T[u-1,,1] <- alphacut(ST, seq(0,1,1/breakpoints))[,1]
T[u-1,,2] <- rev(alphacut(ST, seq(0,1,1/breakpoints))[,2])
H[u-1,,1] <- alphacut(SH, seq(0,1,1/breakpoints))[,1]
H[u-1,,2] <- rev(alphacut(SH, seq(0,1,1/breakpoints))[,2])
df[u-1,] <- r[u-1,] - 1
#E[u-1,] <- c(supp(SE)[1], core(SE)[1], core(SE)[2], supp(SE)[2])
#T[u-1,] <- c(supp(ST)[1], core(ST)[1], core(ST)[2], supp(ST)[2])
#H[u-1,] <- c(supp(SH)[1], core(SH)[1], core(SH)[2], supp(SH)[2])
#df[u-1,] <- r[u-1,] - 1
}
Sum.T <- H+E
#all.equal(T,(H+E))
if(length(grep(':',colnames(data))) != 0){
Ht<- array(rep(0), dim=c(nc,b,2))
SSE <- array(rep(0), dim=c(length(grep(':',colnames(data))),b,2))
SSMAIN <- array(rep(0), dim=c(length(grep(':',colnames(data))),b,2))
Y.Square.obs <- array(rep(0), dim=c(nrow(data),b,2))
for(w in 1:nrow(data)){
#Y.obs <- cbind(Y[w,,1], rev(Y[w,,2])); colnames(Y.obs) <- c("L","U")
#Y.dist[w,1] <- distance(Y.obs, TrapezoidalFuzzyNumber(0,0,0,0), type = distance.type, i=i, j=j, theta = theta, p=p, q=q, breakpoints = breakpoints)
Y.obs <- c(Y[w,1,1], Y[w,(breakpoints+1),1], Y[w,1,2], Y[w,(breakpoints+1),2])
if (is.unsorted(Y.obs) == TRUE){Y.obs <- sort(Y.obs)}
Y.obs <- (PiecewiseLinearFuzzyNumber(Y.obs[1], Y.obs[2], Y.obs[3], Y.obs[4]))*(PiecewiseLinearFuzzyNumber(Y.obs[1], Y.obs[2], Y.obs[3], Y.obs[4]))
Y.Square.obs[w,,1] <- alphacut(Y.obs, seq(0,1,1/breakpoints))[,1]
Y.Square.obs[w,,2] <- rev(alphacut(Y.obs, seq(0,1,1/breakpoints))[,2])
# c(supp(Y.obs)[1], core(Y.obs)[1], core(Y.obs)[2], supp(Y.obs)[2])
}
for(u in grep(':',colnames(data))){
S5 <- 0
S3 <- matrix(rep(0), nrow = b, ncol=2)
S4 <- matrix(rep(0), nrow = b, ncol=2)
SS5 <- matrix(rep(0), nrow = b, ncol=2)
SSInt <- matrix(rep(0), nrow = b, ncol=2)
init <- match(colnames(data)[u], colnames(data)[grep(':',colnames(data))])
j.1 <- match(strsplit(colnames(data)[grep(':',colnames(data))][init], ":")[[1]][1], colnames(data))
k.1 <- match(strsplit(colnames(data)[grep(':',colnames(data))][init], ":")[[1]][2], colnames(data))
S3[,1] <- H[j.1-1,,1] + colSums(Y.Square.obs)[,1]/nt[u-1,1]
S3[,2] <- rev(H[j.1-1,,2] + colSums(Y.Square.obs)[,2]/nt[u-1,1])
S4[,1] <- H[k.1-1,,1] + colSums(Y.Square.obs)[,1]/nt[u-1,1]
S4[,2] <- rev(H[k.1-1,,2] + colSums(Y.Square.obs)[,2]/nt[u-1,1])
for (v in 1:(r[(j.1-1),1]*r[(k.1-1),1])){
mat <- Y.Square.obs[which((data[,j.1]==as.numeric(expand.grid(levels(data[,j.1]), levels(data[,k.1]))[v,1]))&(data[,k.1]==as.numeric(expand.grid(levels(data[,j.1]), levels(data[,k.1]))[v,2]))),,]
SS5[,1] <- colSums(mat)[,1]/nrow(mat)
SS5[,2] <- rev(colSums(mat)[,2]/nrow(mat))
S5 <- S5 + SS5
}
#SSInt <- S5 + sum(Y.dist)^2/nt[u-1,1] - (S3 + S4)
SSquare <- cbind(colSums(Y.Square.obs)[,1], rev(colSums(Y.Square.obs)[,2]))/nt[u-1,1]
SSInt <- Fuzzy.Difference(S5+SSquare,S3+S4, alphacuts = TRUE, breakpoints = breakpoints)
Ht[u-1,,1] <- SSInt[,1]
Ht[u-1,,2] <- rev(SSInt[,2])
#Ht[u-1,1] <- SSInt
#SSInt[,1] <- S5[,1] + colSums(Y.Square.obs)[,1]/nt[u-1,1] - (S3[,1] + S4[,1])
#SSInt[,2] <- S5[,2] + rev(colSums(Y.Square.obs)[,2])/nt[u-1,1] - (S3[,2] + S4[,2])
#SSE[init,1] <- T[u-1,1] + sum(Y.dist)^2/nt[u-1,1] - S5
SSE[init,,1] <- (Fuzzy.Difference(cbind(T[u-1,,1], rev(T[u-1,,2]))+SSquare, S5, alphacuts = TRUE, breakpoints = breakpoints))[,1]
SSE[init,,2] <- rev(Fuzzy.Difference(cbind(T[u-1,,1], rev(T[u-1,,2]))+SSquare, S5, alphacuts = TRUE, breakpoints = breakpoints)[,2])
SSMAIN[init,,1] <- (Fuzzy.Difference(S3+S4, 2*SSquare, alphacuts = TRUE, breakpoints = breakpoints))[,1]
SSMAIN[init,,2] <- rev(Fuzzy.Difference(S3+S4, 2*SSquare, alphacuts = TRUE, breakpoints = breakpoints)[,2])
#SSMAIN[init,1] <- S3 + S4 - 2*(sum(Y.dist)^2/nt[u-1,1])
df[u-1,] <- (r[j.1-1,1] - 1) * (r[k.1-1,1] - 1)
}
H[(grep(':',colnames(data))[1]-1):nc,,1] <- Ht[(grep(':',colnames(data))[1]-1):nc,,1]
H[(grep(':',colnames(data))[1]-1):nc,,2] <- Ht[(grep(':',colnames(data))[1]-1):nc,,2]
}
# Faut attacher les donnees
# pour le cas taille utiliser plutot data <- mf
if (is.balanced(ni[1:(ncol(mf)-1),]) == FALSE){
seq <- SEQ.ORDERING.APPROXIMATION(scope = formula, data = data, f.response = Y)
E[1:(ncol(data)-1),,] <- seq$E.cond
H[1:(ncol(data)-1),,1] <- rbind(H[1,,1],seq$H.cond[,,1])
H[1:(ncol(data)-1),,2] <- rbind(H[1,,2],seq$H.cond[,,2])
coef.model <- coefficients(seq)
predicted.values <- fitted.values(seq)
residuals <- residuals(seq)
} else{
#if ((op.contrasts %in% c("contr.helmert", "contr.poly")) == FALSE) {stop("Error in constrasts")}
coef<- array(rep(0), dim=c(nc,max(r)*b,2))
#contrasts.default <- matrix(rep(0), nrow= nc, ncol = max(r))
for(u in 2:ncol(data)){
#contrasts(data[,u]) <- op.contrasts
# contrasts.default[u-1,1:r[u-1,1]] <- t(contrasts(data[,u]))[1,]
for(v in 1:r[u-1,1]){
Part.mean <- cbind(mat.means[u-1,((b*(v-1)+1):(b*(v-1)+b)),1], rev(mat.means[u-1,((b*(v-1)+1):(b*(v-1)+b)),2]))
colnames(Part.mean) <- c("L", "U")
Y..mean <- cbind(Y..[u-1,,1], rev(Y..[u-1,,2]))
colnames(Y..mean) <- c("L","U")
coef[u-1,((b*(v-1)+1):(b*(v-1)+b)),1] <- Fuzzy.Difference(Part.mean, Y..mean, alphacuts = TRUE, breakpoints = breakpoints)[,1]
coef[u-1,((b*(v-1)+1):(b*(v-1)+b)),2] <- rev(Fuzzy.Difference(Part.mean, Y..mean, alphacuts = TRUE, breakpoints = breakpoints)[,2])
#coef[u-1,v] <- distance(Part.mean, TrapezoidalFuzzyNumber(0,0,0,0), type = distance.type, i=i, j=j, theta = theta, p=p, q=q, breakpoints = breakpoints) -
# distance(TrapezoidalFuzzyNumber(0,0,0,0), Y..mean, type = distance.type, i=i, j=j, theta = theta, p=p, q=q, breakpoints = breakpoints)
}
}
coef.var <- coef
coef.model <- array(rep(0), dim=c((nc + 1),b,2))
coef.model[1,,1] <- Y..[1,,1]
coef.model[1,,2] <- Y..[1,,2]
coef.model[2:(nrow(coef.var)+1),,1] <- coef.var[1:nrow(coef.var),((b+1):(2*b)),1]
coef.model[2:(nrow(coef.var)+1),,2] <- coef.var[1:nrow(coef.var),((b+1):(2*b)),2]
if(length(grep(':',colnames(data))) != 0){
val.int <- t(t(grep(":",colnames(data))))
coef.int <- array(rep(0), dim=c(length(val.int),b,2))
Scoefi1 <- 0
Scoefi2 <- 0
e=1
for(l in val.int[order(-val.int)]){
#coef.int[(length(val.int)-e+1),1] <- sum(coef[l-1,])
for(v in 1:r[u-1,1]){
Scoefi1 <- Scoefi1 + coef[l-1,((b*(v-1)+1):(b*(v-1)+b)),1]
Scoefi2 <- Scoefi2 + coef[l-1,((b*(v-1)+1):(b*(v-1)+b)),2]
}
coef.int[(length(val.int)-e+1),,1] <- Scoefi1
coef.int[(length(val.int)-e+1),,2] <- Scoefi2
coef.var <- coef.var[-(l-1),,]
e=e+1
}
coef.model[(nrow(coef.var)+2):nrow(coef.model),,1] <- coef.int[1:nrow(coef.int),,1]
coef.model[(nrow(coef.var)+2):nrow(coef.model),,2] <- coef.int[1:nrow(coef.int),,2]
}
#coef.model <- c(distance(AY.., TrapezoidalFuzzyNumber(0,0,0,0)),
# coef.var[,2], coef.int)
mat.matrix <- data.frame(as.matrix(model.matrix(formula, data)))
predicted.values <- fuzzy.predicted.values(dataset = mat.matrix, coef.model = coef.model)
#predicted_values <- as.matrix(model.matrix(formula, data)) %*% (coef.model)
residuals <- fuzzy.residuals(data.fuzzified, predicted.values)
}
# Test de Fisher
# --------------
# F-value for the full model
p <- nc
CSH <- cbind(colSums(H[,,1]), rev(colSums(H[,,2])))/p
if (is.alphacuts(CSH)==FALSE){
CSH <- sort(CSH[1,1], CSH[breakpoints+1,1], CSH[1,2], CSH[breakpoints+1,2])
CSH <- TrapezoidalFuzzyNumber(CSH[1], CSH[2], CSH[3], CSH[4])}
CST <- Sum.T[grep(max(Sum.T[,101,2]), Sum.T[,101,2]),,]
CST[,2] <- rev(CST[,2])
ME.full <- (alphacut(Fuzzy.Difference(CST, CSH), seq(0,1,1/breakpoints)))/(max(nt)-1-p)
Ft <- qf(1-sig, df1=r-1, df2=max(nt)-sum(df)-1)#max(nt)-r)
CSE <- ME.full#*Ft
pvalue.manova.model <- pf(CSE[,2], CSH[,2], df1 = p, df2 = max(nt)-sum(df)-1)#max(nt)-1-p)
plot(pvalue.manova.model, seq(0,1,1/breakpoints), 'l')
pvalue.manova <- array(rep(0), dim=c(nrow(H), breakpoints+1, 2))
MH <- H
ME <- E
for(z in 1:nrow(H)){
MH[z,,1] <- H[z,,1] / df[z,1]
MH[z,,2] <- H[z,,2] / df[z,1]
ST <- cbind(Sum.T[z,,1], rev(Sum.T[z,,2]))
if (is.alphacuts(ST)==FALSE){
ST <- sort(ST[1,1], ST[breakpoints+1,1], ST[1,2], ST[breakpoints+1,2])
ST <- TrapezoidalFuzzyNumber(ST[1], ST[2], ST[3], ST[4])}
Ft <- qf(1-sig, df1=p, df2=max(nt)-sum(df)-1)#max(nt)-r)
ME[z,,1] <- (alphacut(Fuzzy.Difference(ST, CSH), seq(0,1,1/breakpoints))[,1]/(max(nt)-1-p))*Ft[1]#Ft[z,1]
ME[z,,2] <- (rev(alphacut(Fuzzy.Difference(ST, CSH), seq(0,1,1/breakpoints))[,2])/(max(nt)-1-p))*Ft[1]#Ft[z,1]
H0 <- cbind(MH[z,,1], rev(MH[z,,2]))
colnames(H0) <- c("L","U")
H1 <- cbind(MH[z,,1]+1, rev(MH[z,,2]+1))
colnames(H1) <- c("L","U")
t <- cbind(ME[z,,1], rev(ME[z,,2]))
colnames(t) <- c("L","U")
t.L <- t[,"L"]
t.U <- t[,"U"]
H0.L <- H0[,"L"]
H0.U <- H0[,"U"]
#if ( t[1,1] >= H0[1,2] ) {
# pvalue.manova[z,,1] = 2*(1-pf( t.U , H0.L , df1 = r-1, df2 = n-r))
#pvalue.manova[z,,2] = 2*(1-pf( sort(t.L, decreasing = TRUE) , sort(H0.U) , df1 = r-1, df2=n-r))
# }
# else if ( t[1,2] <= H0[1,1] ) {
# pvalue.manova[z,,1] = 2* pf( t.U , H0.L , df1 = r-1, df2 = n-r)
#pvalue.manova[z,,2] = 2* pf( sort(t.L, decreasing = TRUE) , sort(H0.U) , df1 = r-1, df2=n-r)
# } else{
#pvalue.manova[z,,2] = pf( sort(t.L, decreasing = TRUE), sort(H0.U), df1 = r-1, df2 =n-r)
# pvalue.manova[z,,1] = pvalue.manova[z,1,2]
# }
#pvalue.manova[z,,1] = (pf( sort(t.L, decreasing = TRUE), sort(H0.L), df1 = r-1, df2 =n-r))
pvalue.manova[z,,1] = (pf( t.L, H0.U, df1 = r-1, df2 =max(nt)-sum(df)-1))#max(nt)-r))
pvalue.manova[z,,2] = rev(pf( sort(t.U, decreasing = TRUE), sort(H0.U), df1 = r-1, df2 =max(nt)-sum(df)-1))#max(nt)-r))max(nt)-r))
}
if (is.na(index.var)==FALSE){
F.MSTR <- MH
F.MSE <- ME
if (plot == TRUE){
plot(F.MSTR[index.var,,1], seq(0,1,1/breakpoints), type='l', xlim=c(min(F.MSTR[index.var,,],F.MSE[index.var,,]), max(F.MSTR[index.var,,],F.MSE[index.var,,])), col = 'blue', xlab = "x", ylab = "alpha", main="Fuzzy decisions - treatments vs. residuals")
#par(new=TRUE)
opar <- par(new=TRUE, no.readonly = TRUE)
on.exit(par(opar))
plot(rev(F.MSTR[index.var,,2]), seq(0,1,1/breakpoints), type='l', xlim=c(min(F.MSTR[index.var,,],F.MSE[index.var,,]), max(F.MSTR[index.var,,],F.MSE[index.var,,])), col = 'blue', xlab = "x", ylab = "alpha")
#par(new=TRUE)
opar <- par(new=TRUE, no.readonly = TRUE)
on.exit(par(opar))
lines(c(F.MSTR[index.var,breakpoints+1,1],F.MSTR[index.var,1,2]), c(1,1), col = "blue")
#par(new=TRUE)
opar <- par(new=TRUE, no.readonly = TRUE)
on.exit(par(opar))
plot(F.MSE[index.var,,1], seq(0,1,1/breakpoints), type='l', xlim=c(min(F.MSTR[index.var,,],F.MSE[index.var,,]), max(F.MSTR[index.var,,],F.MSE[index.var,,])), col = 'red', xlab = "x", ylab = "alpha")
#par(new=TRUE)
opar <- par(new=TRUE, no.readonly = TRUE)
on.exit(par(opar))
plot(rev(F.MSE[index.var,,2]), seq(0,1,1/breakpoints), type='l', xlim=c(min(F.MSTR[index.var,,],F.MSE[index.var,,]), max(F.MSTR[index.var,,],F.MSE[index.var,,])), col = 'red', xlab = "x", ylab = "alpha")
#par(new=TRUE)
opar <- par(new=TRUE, no.readonly = TRUE)
on.exit(par(opar))
lines(c(F.MSE[index.var,breakpoints+1,1],F.MSE[index.var,1,2]), c(1,1), col = "red")
legend("bottomright", legend=c("F.MSE", "F.MSTR"), col=c("red", "blue"), lty=1)
}
# DECISION RULE 1
#################
Surf.MSTR <- abs(integrate.num(cut=F.MSTR[index.var,,1], alpha=seq(0,1,1/breakpoints), int.method)) + abs(integrate.num(cut=F.MSTR[index.var,,2], alpha=seq(0,1,1/breakpoints), int.method))
Surf.MSE <- abs(integrate.num(cut=F.MSE[index.var,,1], alpha=seq(0,1,1/breakpoints), int.method)) + abs(integrate.num(cut=F.MSE[index.var,,2], alpha=seq(0,1,1/breakpoints), int.method))
convicTR <- Surf.MSTR/(Surf.MSE+Surf.MSTR)
convicE <- Surf.MSE/(Surf.MSE+Surf.MSTR)
if(convicTR >= convicE){
decision <- list(noquote(paste0("Variable index: ", index.var, ". Decision: The null hypothesis (H0) is rejected at the ", sig, " significance level. ")),
noquote(paste0(" Degree of conviction (treatments of ",colnames(mf)[2], ") = ", round(convicTR,5), ".")),
noquote(paste0(" Degree of conviction (residuals) ", round(convicE,5), ".")))
} else {
decision <- list(noquote(paste0("Variable index: ", index.var, ". Decision: The null hypothesis (H0) is not rejected at the ", sig, " significance level. ")),
noquote(paste0(" Degree of conviction (treatments of ",colnames(mf)[2], ") = ", round(convicTR,5), ".")),
noquote(paste0(" Degree of conviction (residuals) ", round(convicE,5), ".")))
}
print(decision)
}
resultFMANOVA = list(formula = formula,
terms = colnames(data),
nlevels = r,
rank = nc,
table = ni,
treatments.SSQ = H,
F.coefficients = F,
pvalue.manova = pvalue.manova,
coefficients = coef.model,
pvalues.coefficients = pvalue.manova,
residuals = residuals,
fitted.values = predicted.values,
total.SSQ.model = CST,
df.total = max(nt)-1,
treatments.SSQ.model = CSH,
error.MSSQ.model = CSE,
df.residuals = max(nt)-1-sum(df),
treatment.SSQ.vars = MH,
df.treatments = df,
residuals.SSQ.vars = ME,
F.model <- Ft,
pvalue.model = pvalue.manova.model,
int.res = if(length(grep(':',colnames(data))) != 0){list(error.SSQ.int = SSE, main.SSQ.int = SSMAIN)} else{NULL}
)
}
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.