Nothing
pivotJ = function(M, j, clear=TRUE, eps=1e-8)
{
for (k in j) {
if (any(abs(M[, k]) > eps) > 0) {
Js = which(abs(as.vector(M[, k])) > eps)
pivotRow = M[Js[1], ]/M[Js[1], k]
nJ = length(Js)
if (nJ > 1) for (i in 2:nJ) M[Js[i], ] = M[Js[i] ,] - M[Js[i], k]*pivotRow
if (clear) {
M[Js[1], ] = 0
} else {
M[Js[1], ] = pivotRow
}
}
}
return(M)
}
sortColName = function(ColNames)
{
sCol = strsplit(ColNames, ":")
nc = length(sCol[[1]])
if (nc > 1) {
nr = length(sCol)
tM = matrix(unlist(sCol), nrow=nr, byrow=TRUE)
sN = "order(nchar(tM[,1]), tM[,1]"
for (i in 2:nc) sN = paste0(sN, ", nchar(tM[,",i,"]), tM[,", i,"]")
sN = paste0(sN, ")")
Ord = eval(parse(text=sN))
Res = ColNames[Ord]
} else {
Res = ColNames
}
return(Res)
}
sumANOVA = function(r1, T1, SST, nObs, yName=NULL)
{
if ("(Intercept)" %in% colnames(r1$g2)) {
DF = c(r1$rank - 1, r1$DFr, nObs - 1)
} else {
DF = c(r1$rank, r1$DFr, nObs)
}
SS = c(SST - r1$SSE, r1$SSE, SST)
if (DF[2] > 0) {
MS = c(SS[1:2]/DF[1:2], NA)
} else {
MS = c(SS[1]/DF[1], NA, NA)
}
if (MS[2] > 0 & DF[2] > 0) {
Fval = c(MS[1]/MS[2], NA, NA)
Pval = c(1 - pf(Fval[1], DF[1], DF[2]), NA, NA)
} else {
Fval = rep(NA, 3)
Pval = rep(NA, 3)
}
ANOVA = cbind(DF, SS, MS, Fval, Pval)
colnames(ANOVA) = c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")
if ("(Intercept)" %in% colnames(r1$g2)) {
rownames(ANOVA) = c("MODEL", "RESIDUALS", "CORRECTED TOTAL")
} else {
rownames(ANOVA) = c("MODEL", "RESIDUALS", "UNCORRECTED TOTAL")
}
if (!is.null(T1)) {
rownames(T1) = paste0(" ", rownames(T1))
ANOVA = rbind(ANOVA[1,,drop=FALSE], T1, ANOVA[2:3,])
}
if (!is.null(yName)) {
attr(ANOVA, "heading") = paste("Response :", yName)
}
class(ANOVA) = "anova"
return(ANOVA)
}
sumREG = function(r1, X)
{
np = ncol(X)
Est = r1$coefficients
if (r1$DFr > 0) {
bVar = r1$g2 %*% crossprod(X) %*% t(r1$g2) * r1$SSE/r1$DFr
bSE = sqrt(diag(bVar))
Tval = r1$coefficients/bSE
Pval = 2*(1 - pt(abs(Tval), r1$DFr))
} else {
bSE = NA
Tval = NA
Pval = NA
}
Est[is.na(r1$DFr2)] = NA
bSE[is.na(r1$DFr2)] = NA
Tval[is.na(r1$DFr2) | is.nan(Tval)] = NA
Pval[is.na(r1$DFr2) | is.nan(Pval)] = NA
vESTM = estmb(diag(np), X, r1$g2)
if (sum(vESTM) == np) {
Parameter = cbind(Est, bSE, r1$DFr2, Tval, Pval)
colnames(Parameter) = c("Estimate", "Std. Error", "Df", "t value", "Pr(>|t|)")
} else {
Parameter = cbind(Est, vESTM, bSE, r1$DFr2, Tval, Pval)
colnames(Parameter) = c("Estimate", "Estimable", "Std. Error", "Df", "t value", "Pr(>|t|)")
}
rownames(Parameter) = colnames(X)
class(Parameter) = "anova"
# attr(Parameter, "R2") = r1$R2
# attr(Parameter, "R2ADJ") = r1$R2ADJ
return(Parameter)
}
lsm0 = function(x, rx, Formula, Data, conf.level=0.95, hideNonEst=TRUE)
{
Lx = llsm0(Formula, Data)
b = rx$coefficients[colnames(Lx)]
PE = Lx %*% b
if (rx$DFr > 0) {
Var = Lx %*% rx$g2 %*% t(Lx) * rx$SSE/rx$DFr
SE = sqrt(diag(Var))
} else {
SE = NA
}
DE = qt((1 + conf.level)/2, rx$DFr)*SE
LL = PE - DE
UL = PE + DE
Res = cbind(PE, LL, UL, SE, rx$DFr)
colnames(Res) = c("LSmean", "LowerCL", "UpperCL", "SE", "Df")
if (hideNonEst) Res[attr(Lx, "nonEst"),] = NA
return(Res)
}
llsm0 = function(Formula, Data)
{
# Find continuous variables at model frame
x = ModelMatrix(Formula, Data)
if (!attr(terms(x), "response")) stop("Dependent variable should be provided")
mf = model.frame(Formula, Data)
if (!is.numeric(mf[,1])) stop("Dependent variable should be numeric!")
cn0 = colnames(mf)
nc0 = ncol(mf)
vc0 = rep(FALSE, nc0) # vector continuous 0
for (i in 1:nc0) vc0[i] = is.numeric(mf[,i])
cv0 = cn0[vc0][-1] # continuous x variable names
ncv0 = length(cv0)
Labels = labels(terms(x))
nLabel = length(Labels)
tN = attr(terms(x), "factors") # table for nest information
fIntercept = attr(x$terms, "intercept") # intercept of original formula
if (fIntercept) { # if original formula does not have intercept, add intercept
X = x$X
} else {
X = cbind(1, x$X)
colnames(X) = c("(Intercept)", colnames(x$X))
x$X = X
x$assign = c(0, x$assign)
for (i in 1:nLabel) x$termIndices[[i]] = x$termIndices[[i]] + 1
attr(x$terms, "intercept") = 1
}
XpX = crossprod(X)
nc = ncol(XpX)
L0 = matrix(NA, nrow=nc, ncol=nc)
dimnames(L0) = dimnames(XpX)
L0[,1] = 1 # intercept column
L1 = as.numeric(XpX[1,] > 0) # intercept row, columns with observation count > 0
# First row (intercept)
for (i in 1:nLabel) {
cI = x$termIndices[Labels[i]][[1]]
if (sum(tN[,Labels[i]] > 1) > 0) { # nested
nestF = rownames(tN)[tN[,Labels[i]] > 1]
nestF = paste(nestF, collapse=":")
nfCI = x$termIndices[[nestF]]
nnfCI = length(nfCI)
nCI = length(cI)
sG = nCI/nnfCI
if (sG > 1) {
for (j in 1:nnfCI) {
ccI = cI[((j - 1)*sG + 1):(j*sG)]
L0[1, ccI] = 1/sum(L1[ccI])*L0[1, nfCI[j]]
}
} else {
L0[1, cI] = 1/sum(L1[cI])
}
} else {
L0[1, cI] = 1/sum(L1[cI])
}
}
for (i in 1:nLabel) {
for (j in 1:nLabel) {
Ls = bL(i, j, x, L0[1,,drop=FALSE])
L0[rownames(Ls), colnames(Ls)] = Ls
}
}
# Fill single continuous columns with mean
if (ncv0 > 0) {
for (i in 1:ncv0) L0[, cv0[i]] = mean(mf[, cv0[i]])
}
# Two or more interaction columns with continuous variable
cn1 = colnames(XpX)
scn1 = strsplit(cn1, ":")
names(scn1) = cn1
for (i in 1:nc) {
ccn = scn1[[i]] # splitted current column name
nccn = length(ccn) # number of ccn
if (nccn == 1) next # already treated above, no more need to consider for single name
cap = intersect(cv0, ccn)
ncap = length(cap)
if (ncap == 0) next # no intersection -> no contiuous varaible
tv = L0[, ccn[1]] # temporary vector
for (j in 2:nccn) tv = tv*L0[, ccn[j]] # it must be more than 1 column (nccn >= 2)
L0[, i] = tv
}
# Check estimability
iNE = (diag(XpX) == 0) # index of not estimable
nNE = sum(iNE)
if (nNE > 0) {
sNE = strsplit(cn1[iNE], ":")
for (i in 1:nNE) {
cStr = sNE[[i]]
for (j in 2:nc) {
cCol = scn1[[j]]
# if (all(cCol %in% cStr) & length(cCol) < length(cStr) & tN[x$assign[j] + 1, x$assign[iNE][i]] == 1) iNE[j] = TRUE
C1 = all(cCol %in% cStr) & length(cCol) < length(cStr)
C2 = length(cCol)
if (C2 > 1) {
if (C1) iNE[j] = TRUE
} else if (C2 == 1) {
if (C1 & tN[x$assign[j] + 1, x$assign[iNE][i]] < 2) iNE[j] = TRUE
}
}
}
}
L0[, diag(XpX) == 0] = 0 # No observation columns
attr(L0, "nonEst") = iNE
if (!fIntercept) { # if original formula does not have intercept
L0 = L0[-1, -1]
attr(L0, "nonEst") = iNE[-1]
}
return(L0)
}
bL = function(m, n, x, L1)
{
Labels = labels(terms(x))
tLabel = Labels[m]
cLabel = Labels[n]
sLabel = strsplit(Labels, ":")
fCap = intersect(sLabel[[m]], sLabel[[n]])
XpX = crossprod(x$X)
cn = colnames(XpX)
iNO = (1:ncol(XpX))[diag(XpX) == 0]
cI = x$termIndices[[tLabel]]
cJ = setdiff(x$termIndices[[cLabel]], iNO)
nr = length(cI)
nc = length(cJ)
Lx = matrix(0, nrow=nr, ncol=nc)
dimnames(Lx) = list(cn[cI], cn[cJ])
if (length(fCap) == 0) {
for (i in 1:nr) Lx[i, ] = L1[1, setdiff(cJ, iNO)]
} else if (m == n) {
Lx = diag(1, nrow=length(cI))
dimnames(Lx) = list(cn[cI], cn[cI])
} else {
tLevels = strsplit(cn[cI], ":")
cLevels = strsplit(cn[cJ], ":")
iLevels = intersect(unlist(tLevels), unlist(cLevels))
tNames = rep("", nr)
cNames = rep("", nc)
for (i in 1:nr) tNames[i] = paste(intersect(tLevels[[i]], iLevels), collapse=":")
for (j in 1:nc) cNames[j] = paste(intersect(cLevels[[j]], iLevels), collapse=":")
for (i in 1:nr) for (j in 1:nc) if (tNames[i] == cNames[j]) Lx[i, j] = 1
Lx = Lx/rowSums(Lx)
}
return(Lx)
}
GrpCode = function(NamesInOrder, rPDIFF, conf.level=0.95)
{
nLevel = length(NamesInOrder)
nL = nrow(rPDIFF)
rowNames = rownames(rPDIFF)
Names = strsplit(rowNames, " - ")
m1 = matrix(NA, nrow=nLevel, ncol=nLevel)
colnames(m1) = NamesInOrder
rownames(m1) = NamesInOrder
iPr = ncol(rPDIFF)
for (k in 1:nL) {
cX = Names[[k]][1]
cY = Names[[k]][2]
if (rPDIFF[k, iPr] < (1 - conf.level)) {
m1[cX, cY] = 0
m1[cY, cX] = 0
}
}
cI = 1
for (j in 1:nLevel) {
fUsed = FALSE
for (i in j:nLevel) {
if (is.na(m1[i, j])) {
m1[i,j] = cI
fUsed = TRUE
} else {
break
}
}
if (j > 1) {
if (m1[i, j - 1] > 0) m1[j:nLevel, j] = m1[i, j - 1]
}
if (fUsed) cI = cI + 1
}
m1[is.na(m1)] = 0
for (j in 2:nLevel) {
fInc = TRUE
for (i in j:nLevel) {
if (m1[i,j] > 0 & m1[i, j - 1] == 0) fInc = FALSE
}
if (fInc) {
for (i in j:nLevel) if (m1[i,j] > 0) m1[i,j] = m1[i, j - 1]
} else {
for (i in j:nLevel) if (m1[i,j] > 0) m1[i,j] = m1[j - 1, j - 1] + 1
}
}
gCode = rep("", nLevel)
for (i in 1:nLevel) {
tv = sort(unique(m1[i,]))
t1 = paste(rep(" ", min(tv[tv > 0]) - 1), collapse="")
t2 = paste(unique(LETTERS[m1[i,1:i]]), collapse="")
t3 = paste(rep(" ", max(m1) - nchar(t1) - nchar(t2)), collapse="")
gCode[i] = paste0(t1, t2, t3)
}
return(gCode)
}
plotDiff = function(lSM, m0, conf.level=0.95, ...)
{
nL = nrow(m0)
nc = ncol(m0)
nLevel = length(lSM)
rowNames = rownames(m0)
Names = strsplit(rowNames, " - ")
Lnames = sort(unique(unlist(Names)))
hDL0 = max((m0[,3] - m0[,2])/4) # Half DL is necessary
if (is.na(hDL0) | is.nan(hDL0)) stop("SE is not available!")
xmin = min(lSM) - hDL0
xmax = max(lSM) + hDL0
Args = list(...)
nArgs = names(Args)
Args$x = 0
Args$y = 0
Args$xlim = c(xmin, xmax)
Args$ylim = c(xmin, xmax)
Args$type = "n"
if ("ref" %in% nArgs) Args$ref = NULL # remove ref argument before calling plot
if (!("xlab" %in% nArgs)) Args$xlab = "Levels"
if (!("ylab" %in% nArgs)) Args$ylab = "Difference from the reference level"
do.call(plot, Args)
abline(a=0, b=1, lty=2)
abline(h=lSM, lty=3)
abline(v=lSM, lty=3)
text(x=xmax, y=lSM, labels=names(lSM))
text(x=lSM, y=xmin, labels=names(lSM))
m2Col = c("x", "y", "hDL")
m2 = matrix(nrow=nL, ncol=length(m2Col))
colnames(m2) = m2Col
for (i in 1:nL) {
m2[i,"x"] = lSM[Names[[i]][1]][[1]]
m2[i,"y"] = lSM[Names[[i]][2]][[1]]
}
m2[,"hDL"] = (m0[,3] - m0[,2])/4 # PE -> PE/sqrt(2) -> DL/PE/sqrt(2) -> DL/PE/sqrt(2)/sqrt(2)*PE -> DL/2
for (i in 1:nL) {
if (m0[i, 1] > 0) { # For left upper region plot, PE of difference should be negative
temp = m2[i, "x"]
m2[i, "x"] = m2[i, "y"]
m2[i, "y"] = temp
}
hDL = m2[i,"hDL"] # Half DL, -> DL/sqrt(2) if 45 degree rotation
if (m0[i, nc] < 1 - conf.level) { Col="blue"
} else { Col="red" }
points(m2[i,"x"], m2[i,"y"], pch=16, col=Col)
lines(c(m2[i,"x"] - hDL, m2[i,"x"] + hDL), c(m2[i,"y"] + hDL, m2[i,"y"] - hDL), col=Col, lwd=2)
}
}
plotDunnett = function(m0, ...)
{
nL = nrow(m0)
rowNames = rownames(m0)
Names = strsplit(rowNames, " - ")
xNames = vector(length=nL)
for (i in 1:nL) xNames[i] = Names[[i]][1]
ymax = max(abs(m0[,2:3]))
Args = list(...)
nArgs = names(Args)
Args$x = 0
Args$y = 0
Args$xlim = c(0.5, nL + 0.5)
Args$ylim = c(-ymax, ymax)
Args$type = "n"
Args$xaxt = "n"
if ("ref" %in% nArgs) Args$ref = NULL # remove ref argument before calling plot
if (!("xlab" %in% nArgs)) Args$xlab = "Levels"
if (!("ylab" %in% nArgs)) Args$ylab = "Difference from the reference level"
if (!("main" %in% nArgs)) Args$main = paste("Difference from the reference level:", Names[[1]][2])
do.call(plot, Args)
abline(h=0)
axis(1, at=1:nL, labels=xNames)
points(x=1:nL, y=m0[,1], pch=16)
for (i in 1:nL) arrows(x0=i, y0=m0[i,2], x1=i, y1=m0[i,3], length=0.1, angle=90, code=3)
}
CheckAlias = function(Formula, Data) # Why did I use this instead of 'alias' function? -> Maybe alias() requires MASS package.
{ # 'alias' seems more pecise, because it uses fractional number (i.e. a/b).
Res = list(Model = Formula)
QR = qr(model.matrix(Formula, Data))
Rank = QR$rank
np = dim(QR$qr)[2]
if (is.null(np) || Rank == np) {
Aliased = NULL
} else {
ip = 1:Rank
dn = dimnames(QR$qr)[[2]]
Aliased = backsolve(QR$qr[ip, ip], QR$qr[ip, -ip, drop=F])
dimnames(Aliased) = list(dn[ip], dn[-ip])
Aliased = t(zapsmall(Aliased))
}
Res$Complete = Aliased
return(Res)
}
ex = function(x1, x2, g2, eps=1e-8)
{
eps2 = eps/max(abs(x1$X))
Res = NULL
L1 = crossprod(x1$X)
nc = NCOL(L1)
if (nc > 1) {
for (i in 1:(nc - 1)) {
if (abs(L1[i, i]) < eps) { L1[i, ] = 0 ; L1[i:nc, i] = 0 ; next }
for (j in (i + 1):nc) L1[j, ] = L1[j, ] - L1[j, i]/L1[i, i]*L1[i, ]
if (abs(L1[i, i]) > eps) L1[i, ] = L1[i, ]/L1[i, i]
}
} else {
L1[1, 1] = 1
}
rownames(L1) = paste0("L", 1:nrow(L1))
L1[abs(L1) < eps2] = 0
Res$e1 = L1
XpX = crossprod(x2$X)
M0 = G2SWEEP(XpX, Augmented=FALSE, eps=eps) %*% XpX
M0[abs(M0) < eps2] = 0
rownames(M0) = paste0("L", 1:nc)
Labels = labels(terms(x2))
nLabel = length(Labels)
LLabel = strsplit(Labels, ":")
if (attr(x2$terms, "intercept")) { re2 = c(1, rep(0, nc - 1))
} else { re2 = NULL }
re3 = NULL
for (i in 1:nLabel) {
Label1 = Labels[i]
Label2 = NULL
for (j in 1:nLabel) {
if (i != j & all(LLabel[[i]] %in% LLabel[[j]])) Label2 = c(Label2, Labels[j])
}
Col1 = x2$termIndices[Label1][[1]]
Col2 = NULL
if (length(Label2) > 0) {
for (j in 1:length(Label2)) Col2 = c(Col2, x2$termIndices[Label2[j]][[1]])
}
Col0 = setdiff(1:nc, c(Col1, Col2))
X0 = x2$X[, Col0]
X1 = x2$X[, Col1]
X2 = x2$X[, Col2]
if (NCOL(X0) > 0) {
Mx = X0 %*% G2SWEEP(crossprod(X0), Augmented=FALSE, eps=eps) %*% t(X0)
Mx[abs(Mx) < eps] = 0
M = diag(NCOL(Mx)) - Mx
X1pM = crossprod(X1, M)
} else {
X1pM = t(X1)
}
X1pMX1 = X1pM %*% X1
gX1pMX1 = G2SWEEP(X1pMX1, Augmented=FALSE, eps=eps)
L2 = M0[x2$termIndices[Label1][[1]], , drop=FALSE]
L2[, Col0] = 0
L2[, Col1] = gX1pMX1 %*% X1pMX1
L2[, Col2] = gX1pMX1 %*% X1pM %*% X2
L2[abs(L2) < eps] = 0
DF1 = NROW(L2[!apply(L2, 1, function(x) all(abs(x) < eps)), , drop=FALSE])
R1 = t(M0[Col1, , drop=FALSE])
R2 = t(M0[Col2, , drop=FALSE])
L3 = t(R1 - R2 %*% G2SWEEP(crossprod(R2)) %*% crossprod(R2, R1))
if (NROW(L3) > 0) {
L3 = pivotJ(L3, Col0, clear=TRUE)
L3 = L3[!apply(L3, 1, function(x) all(abs(x) < eps)), , drop=FALSE]
}
DF2 = NROW(L3)
if (DF2 != DF1) {
if (!is.null(Col2)) {
Ms = pivotJ(M0[c(Col1, Col2), , drop=FALSE], Col0, clear=TRUE)
Ms = pivotJ(Ms, Col1, clear=FALSE)
fbazr = apply(abs(Ms[, Col1, drop=FALSE]), 1, max) < eps
bazr = which(fbazr)
bnazr = which(!fbazr)
if (length(bnazr) == 0) {
L3 = NULL
} else if (length(bnazr) <= length(Col1)) {
R1 = t(Ms[bnazr, , drop=FALSE])
R2 = t(Ms[bazr, , drop=FALSE])
L3 = t(R1 - R2 %*% G2SWEEP(crossprod(R2)) %*% crossprod(R2, R1))
rownames(L3) = (rownames(M0)[Col1])[1:NROW(L3)]
} else {
R1 = t(M0[Col1, , drop=FALSE])
R2 = t(M0[Col2, , drop=FALSE])
L3 = t(R1 - R2 %*% G2SWEEP(crossprod(R2)) %*% crossprod(R2, R1))
}
} else {
L3 = pivotJ(M0[Col1, , drop=FALSE], Col1, clear=FALSE)
L3 = pivotJ(L3, Col0, clear=TRUE)
}
if (NROW(L3) > 0) {
L3 = L3[!apply(L3, 1, function(x) all(abs(x) < eps)), , drop=FALSE]
}
DF3 = NROW(L3)
} else {
DF3 = DF1
}
if (DF3 == 0) L2 = matrix(0, nrow(L2), ncol(L2))
re2 = rbind(re2, L2)
if (!is.null(L3)) re3 = rbind(re3, L3)
}
rownames(re2) = paste0("L", 1:NCOL(re2))
re2[abs(re2) < eps2] = 0
Res$e2 = re2
M0[, ] = 0 # clear all
if (attr(x2$terms, "intercept")) M0[1, 1] = 1 # Intercept
if (!is.null(re3)) M0[rownames(re3), ] = re3
M0[abs(M0) < eps2] = 0
Res$e3 = M0
return(Res) # Do not use zapsmall !
}
WhiteH = function(X, we2) # we2 = weight * e^2
{
n = NROW(X)
K = NCOL(X)
Df = K*(K + 1)/2
e2.in = scale(we2, scale=F)
psi.i = matrix(nrow=n, ncol=Df)
for (i in 1:n) psi.i[i, ] = crossprod(X[i, , drop=F])[lower.tri(diag(K), T)]
psi.in = scale(psi.i, scale=F)
Dn = apply(psi.in, 2, function(x) mean(e2.in*x))
sumB = matrix(0, nrow=Df, ncol=Df)
for (i in 1:n) sumB = sumB + e2.in[i]^2*crossprod(psi.in[i, , drop=F])
iBn = G2SWEEP(sumB/n)
Df1 = attr(iBn, "rank")
White = n*t(Dn) %*% iBn %*% Dn
p.val = 1 - pchisq(White, Df1)
return(c(Chisq=White, Df=Df1, p=p.val))
}
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.