set.seed(1)
## Scale yields to reproduce analyses reported in Vargas et al (2001)
yield.scaled <- wheat$yield * sqrt(3/1000)
## Reproduce (up to error caused by rounding) Table 1 of Vargas et al (2001)
aov(yield.scaled ~ year*tillage*summerCrop*manure*N, data = wheat)
treatment <- interaction(wheat$tillage, wheat$summerCrop, wheat$manure,
wheat$N, sep = "")
mainEffects <- lm(yield.scaled ~ year + treatment, data = wheat)
svdStart <- residSVD(mainEffects, year, treatment, 3)
bilinear1 <- update(asGnm(mainEffects), . ~ . +
Mult(year, treatment),
start = c(coef(mainEffects), svdStart[,1]))
bilinear2 <- update(bilinear1, . ~ . +
Mult(year, treatment, inst = 2),
start = c(coef(bilinear1), svdStart[,2]))
bilinear3 <- update(bilinear2, . ~ . +
Mult(year, treatment, inst = 3),
start = c(coef(bilinear2), svdStart[,3]))
anova(mainEffects, bilinear1, bilinear2, bilinear3)
## Examine the extent to which, say, mTF explains the first bilinear term
bilinear1mTF <- gnm(yield.scaled ~ year + treatment + Mult(1 + mTF, treatment),
family = gaussian, data = wheat)
anova(mainEffects, bilinear1mTF, bilinear1)
## How to get the standard SVD representation of an AMMI-n model
##
## We'll work with the AMMI-2 model, which here is called "bilinear2"
##
## First, extract the contributions of the 5 terms in the model:
##
wheat.terms <- termPredictors(bilinear2)
##
## That's a matrix, whose 4th and 5th columns are the interaction terms
##
## Combine those two interaction terms, to get the total estimated
## interaction effect:
##
wheat.interaction <- wheat.terms[, 4] + wheat.terms[, 5]
##
## That's a vector, so we need to re-shape it as a 24 by 10 matrix
## ready for calculating the SVD:
##
wheat.interaction <- matrix(wheat.interaction, 24, 10)
##
## Now we can compute the SVD:
##
wheat.interaction.SVD <- svd(wheat.interaction)
##
## Only the first two singular values are nonzero, as expected
## (since this is an AMMI-2 model, the interaction has rank 2)
##
## So the result object can be simplified by re-calculating the SVD with
## just two dimensions:
##
wheat.interaction.SVD <- svd(wheat.interaction, nu = 2, nv = 2)
Run the code above in your browser using DataLab