# R code to fit an Aitken model
# the metabolite data from section 14 of the notes
genotype <-c('A','B','B','A','A','B','B','A')
gen <- as.factor(genotype)
nsdl <- c(9,6,7,9,6,6,8,5)
metab <- c(0.0231,0.0275,0.0521,0.0001,0.0063,0.0138,0.1061,0.0482)
# the ols analysis
ols.lm <- lm(metab~gen)
summary(ols.lm)
# an Aitken model analysis, by hand
# v is an 8 x 8 matrix with 1/nsdl on the diagonal
v <- diag(1/nsdl)
# need to get the inverse square root matrix of v
# use eigen function to do that
temp <- eigen(v)
# returns a list with two components: values, a vector, and vectors, a matrix
u <- temp$vectors
d <- diag(temp$values) # convert e-vals to a diagonal matrix
round(v - u %*% d %*% t(u), 5) # check u d t' = v
all(v == u %*% d %*% t(u)) # another possible check, more sensitive to num. error
svi <- u %*% diag(1/sqrt(temp$values)) %*% t(u) # inverse square root matrix
# check that works the way we want it to
v %*% svi %*% t(svi) # svi svi' = v^{-1}, and v * v^{-1} = I
z <- svi %*% metab
w <- svi %*% model.matrix(ols.lm) # could also use any other equivalent X matrix
gls.lm <- lm(z~-1+w) # need to suppress default intercept
summary(gls.lm)
# could also do using weighted least squares, since V diagonal
# weights are 1/variance
wls.lm <- lm(metab~gen, weight=nsdl)
summary(wls.lm)
# some useful hints / tricks to construct different sorts of V matrices
# AR(1) structure. V has bands. This sort of matrix is a toeplitz matrix
toeplitz(1:5)
# to general a 8 x 8 ar(1) matrix with specified rho
rho <- 0.7
v <- toeplitz(rho^(0:7)) # sequence starts at 0 so diag = 1, ends at rho^7
# compound symmetry structure (obs within groups), e.g. example 4
grp <- c(1,1,2,2,3,3,4,4,4) # or use rep(1:4,c(2,2,2,3))
v <- outer(grp, grp, '==') # true (i.e. 1) if grp[i] = grp[j]
k <- 2/1.5 # ratio of sigma^2_p / sigma^2_c
v <- v + diag(rep(k,length(grp))) # add k to diagonal, has side effect of coercing T/F to 1/0