# SSE via squared distances (for procD.lm RRPP method et al.) SSE <- function(L){# L is a linear model r <- as.matrix(resid(L)) S <- r%*%t(r) sse <- sum(diag(S)) sse } # P-values (for procD.lm RRPP method) pval = function(s){# s = sampling distribution p = length(s) r = rank(s)[1]-1 pv = 1-r/p pv } # ProcD.lm modified to do analysis through the origin procD.pic <- function(f1, iter = 999){ data=NULL form.in <- formula(f1) Terms <- terms(form.in, keep.order = TRUE) Y <- eval(form.in[[2]], parent.frame()) df <- df.tmp <- SS.tmp <- SS.obs <- F <- Rsq <- array() dat <- model.frame(form.in, data=NULL) Xdims <- dim(as.matrix(model.matrix(Terms))) j <- ncol(attr(Terms, "factors")) if(j==1){ mod.mat <- model.matrix(Terms[1], data = dat) SS.res <- SSE(lm(Y~mod.mat[,-1]-1)) #NOTE, This is correct RESID SS! predmod<-predict(lm(Y~mod.mat[,-1]-1)) #brute-force SSModel. Not elegant, but works. SS.obs<-sum(diag(predmod%*%t(predmod))) SS.tot <- SS.obs+SS.res df <- ncol(mod.mat)-1 MS <- SS.obs/df df.tot <- nrow(Y) #NOTE: not n-1 b/c PIC are already N-1. Otherwise, F ratios won't perfectly match df.res <- nrow(Y) - df ## Ditto here, forget the -1 MS.tot <- SS.tot/df.tot MS.res <- SS.res/df.res Rsq <- SS.obs/SS.tot F <- MS/MS.res P <- array(c(SS.obs,rep(0,iter))) for(i in 1:iter){ Yr <- Y[sample(nrow(Y)),] predmod.r<-predict(lm(Yr~mod.mat[,-1]-1)) #brute-force SSModel. Not elegant, but works. P[i+1]<-sum(diag(predmod.r%*%t(predmod.r))) # P[i+1] <- SS.null-SS.tmp } P.val <- pval(P) } ##NOTE: haven't fixed factorial version yet.... if(j>1){ Xs <- array(0, c(Xdims, (j+1))) Xs[,1,1] <- 1 for(i in 1:j){ x <- as.matrix(model.matrix(Terms[1:i], data=dat)) Xs[,1:ncol(x),(i+1)] <- x SS.tmp[i] <- SSE(lm(Y ~ x -1)) df.tmp[i] <- ifelse(ncol(x) == 1, 1, (ncol(x) - 1)) ifelse(i == 1, df[i] <- df.tmp[i], df[i] <- (df.tmp[i] - df.tmp[i - 1])) } SS.null <- (c(SSE(lm(Y~1)),SS.tmp))[1:j] SS.obs <- SS.null - SS.tmp MS <- SS.obs/df SS.tot <- SSE(lm(Y~1)) SS.res <- SS.tot - sum(SS.obs) df.tot <- nrow(Y) - 1 df.res <- nrow(Y) - 1 - sum(df) MS.tot <- SS.tot/df.tot MS.res <- SS.res/df.res Rsq <- SS.obs/SS.tot F <- MS/MS.res P <- array(0,c(dim(matrix(SS.obs)),iter+1)) P[,,1] <- SS.obs for(i in 1:iter){ Yr <- Y[sample(nrow(Y)),] for (ii in 1:j) { SS.tmp[ii] <- SSE(lm(Yr ~ Xs[,,ii+1] -1)) SS.null <- (c(SSE(lm(Y~1)),SS.tmp))[1:j] } P[,,i+1] <- SS.null-SS.tmp } P.val <- Pval.matrix(P) } anova.tab <- data.frame(df = c(df,df.res,df.tot), SS = c(SS.obs, SS.res, SS.tot), MS = c(MS, MS.res, MS.tot), Rsq = c(Rsq, NA, NA), F = c(F, NA, NA), P.val = c(P.val, NA, NA)) rownames(anova.tab) <- c(attr(Terms, "term.labels"), "Residuals","Total") anova.title = "\nRandomization of Raw Values used\n" attr(anova.tab, "heading") <- paste("\nType I (Sequential) Sums of Squares and Cross-products\n",anova.title) class(anova.tab) <- c("anova", class(anova.tab)) return(anova.tab) }