## Read data fake.data <- read.table(file="fake.txt", header = T) ## Number of treatments t <- length(levels(fake.data\$Treatment)) ## Number of samples N <- nrow(fake.data) ## The log-counts and 0 for last equation y <- as.vector(t(fake.data["Response"])) ## Assign contrasts to a factor contrasts(fake.data\$Treatment) <- contr.treatment(t) #contrasts(fake.data\$Treatment) <- contr.helmert(t) ## We construct the model ## y[i,j] = mu[i] + e[i,j] ## side condition: ## sum(mu[i]) = 0 ## Full design matrix full.dm <- model.matrix(Response ~ Treatment, fake.data) ## Reduced design matrix reduced.dm <- model.matrix(Response ~ Treatment - Treatment, fake.data) ## Now we fit the data ## ## y = dm %*% x + e, ## ## where x is the vector c(mu, mu[1],...) ## and e_|_ ColSpace(dm). ## ## If dm = Q %*% R then we multiply on the left by t(Q): ## ## R %*% x = t(Q) %*% y ## ## which is triangular. design.fit <- function(design.matrix, y) { qr <- qr(design.matrix) Q <- qr.Q(qr) R <- qr.R(qr) Ri <- solve(R) # Find inverse ## Environment to put variables in env <- new.env() env\$coef <- Ri %*% crossprod(Q, y) # coefficients env\$err <- qr.resid(qr, y) # Residues env\$resid <- sum(env\$err^2) return(as.list(env)) } full <- design.fit(full.dm, y) reduced <- design.fit(reduced.dm, y) SSTotal <- sum((y-mean(y))^2) SSMean <- SSTotal / (N-1) #same as var(y) SSE <- full\$resid MSE <- SSE / (N-t) SST <- reduced\$resid - full\$resid MST <- SST / (t-1) F <- MST / MSE fake.aov <- aov(Response ~ Treatment, fake.data, projections = T)