> <\body> |||<\author-address> Copyright Marek Rychlik, 2009 >> The measured variate is grain yield, subject to several treatment levels of a fertilizer. In this experimental study, the treatment is a quantitative factor (fertilizer level). This is an factor. R treats ordered factors differently from unordered factors: it automatically assigns contrasts, which are discrete orthogonal polynomials. Manually, these contrasts can be generated by calling . <\input| >> t = 5; r = 3; N = r * t <\input| >> Yield = as.vector(matrix(c(12.2, 16.0, 18.6, 17.6, 18.0, \ \ 11.4, 15.5, 20.2, 19.3, 16.4, \ \ 12.4, 16.5, 18.2, 17.1, 16.6), byrow = T, nrow = r, ncol = t)) <\input| >> PlantDensity = gl(t, r, labels = seq(10, 50, 10), ordered = T) <\input| >> PlantDensityLevels = as.numeric(levels(PlantDensity)) <\input| >> yield.data = data.frame(PlantDensity, Yield) <\input| >> yield.data <\output> \ \ \ PlantDensity Yield 1 \ \ \ \ \ \ \ \ \ \ \ 10 \ 12.2 2 \ \ \ \ \ \ \ \ \ \ \ 10 \ 11.4 3 \ \ \ \ \ \ \ \ \ \ \ 10 \ 12.4 4 \ \ \ \ \ \ \ \ \ \ \ 20 \ 16.0 5 \ \ \ \ \ \ \ \ \ \ \ 20 \ 15.5 6 \ \ \ \ \ \ \ \ \ \ \ 20 \ 16.5 7 \ \ \ \ \ \ \ \ \ \ \ 30 \ 18.6 8 \ \ \ \ \ \ \ \ \ \ \ 30 \ 20.2 9 \ \ \ \ \ \ \ \ \ \ \ 30 \ 18.2 10 \ \ \ \ \ \ \ \ \ \ 40 \ 17.6 11 \ \ \ \ \ \ \ \ \ \ 40 \ 19.3 12 \ \ \ \ \ \ \ \ \ \ 40 \ 17.1 13 \ \ \ \ \ \ \ \ \ \ 50 \ 18.0 14 \ \ \ \ \ \ \ \ \ \ 50 \ 16.4 15 \ \ \ \ \ \ \ \ \ \ 50 \ 16.6 <\input| >> \; > plotting parameters> <\input| >> opts = options(); opts\$texmacs\$width=7.5; opts\$texmacs\$height=7.5; opts\$texmacs\$nox11=F; options(opts) <\input| >> \; > <\input| >> X11(pointsize=8, height=3, width=3) <\input| >> \; > <\input| >> xadd = 1; yadd = 1 <\input| >> plot(numeric(0), numeric(0), xlim = range(PlantDensityLevels)+c(-1,1)*xadd, ylim = range(Yield)+c(-1,1)*yadd, type = "n", ann = FALSE) <\input| >> points(PlantDensityLevels[as.numeric(PlantDensity)], Yield); v() <\output> |ps>||||||> <\input| >> \; > This is standard: <\input| >> yield.aov = aov(Yield ~ PlantDensity, yield.data, projections = T, qr = T) <\input| >> summary(yield.aov) <\output> \ \ \ \ \ \ \ \ \ \ \ \ \ Df Sum Sq Mean Sq F value \ \ \ Pr(\F) \ \ \ PlantDensity \ 4 87.600 \ 21.900 \ 29.278 1.690e-05 *** Residuals \ \ \ 10 \ 7.480 \ \ 0.748 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \; --- Signif. codes: \ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\ <\input| >> \; > In this discussion, we utilize the function which does the heavy lifting of decomposing experimental error and between-groups variation along the contrast directions. sums of squares used in our analysis are the sums of squares of the columns of the matrix component of the data structure returned by . <\input| >> yield.proj = proj(yield.aov, onedf = T) <\input| >> SS = as.matrix(apply(yield.proj^2, 2, sum)) <\input| >> (SSC = SS[2:t,]) <\output> PlantDensity.L PlantDensity.Q PlantDensity.C PlantDensity^4\ \ \ \ \ \ \ \ \ \ \ 43.2 \ \ \ \ \ \ \ \ \ \ 42.0 \ \ \ \ \ \ \ \ \ \ \ 0.3 \ \ \ \ \ \ \ \ \ \ \ 2.1\ <\input| >> (SST = sum(SSC)) <\output> [1] 87.6 <\input| >> (SSE = SS[t+1,]) <\output> Residuals\ \ \ \ \ \ 7.48\ <\input| >> \; > We utilize the explicit expressions for and other sums of squares. We process all contrasts in parallel, taking advantage of processing of vector components in parallel, offered by the R arithmetic operator and functions. <\input| >> MSC = SSC <\input| >> MSE = SSE / (N - t) <\input| >> (F.contrasts = MSC / MSE) <\output> PlantDensity.L PlantDensity.Q PlantDensity.C PlantDensity^4\ \ \ \ \ 57.7540107 \ \ \ \ 56.1497326 \ \ \ \ \ 0.4010695 \ \ \ \ \ 2.8074866\ <\input| >> MST = SST / (t - 1)\ <\input| >> (F.omnibus = MST / MSE) <\output> Residuals\ \ 29.27807\ <\input| >> \; > We obtain all significance level simultaneously. <\input| >> (sig.contrasts = pf(F.contrasts, df1 = 1, df2 = N - t, lower.tail = F)) <\output> PlantDensity.L PlantDensity.Q PlantDensity.C PlantDensity^4\ \ \ 1.840830e-05 \ \ 2.078706e-05 \ \ 5.407497e-01 \ \ 1.247622e-01\ <\input| >> (sig.omnibus = pf(F.omnibus, df1 = t - 1, df2 = N - t, lower.tail = F)) <\output> \ \ \ Residuals\ 1.689528e-05\ <\input| >> \; > We can see that our results, within numerical precision, agree with the results of the omnibus F-statistic summarized by and . <\input| >> summary(yield.aov) <\output> \ \ \ \ \ \ \ \ \ \ \ \ \ Df Sum Sq Mean Sq F value \ \ \ Pr(\F) \ \ \ PlantDensity \ 4 87.600 \ 21.900 \ 29.278 1.690e-05 *** Residuals \ \ \ 10 \ 7.480 \ \ 0.748 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \; --- Signif. codes: \ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\ <\input| >> \; > The coefficients of the polynomial approximation are obtained by calling on the result of . The functions and are used for the following two tasks: <\itemize> To obtain the values of the orthogonal polynomials at the abscissas which are the levels of the factor To interpolate the values of the orthogonal polynomials for a much finer mesh, on which we would like to approximate the approximating polynomial The evaluation of a polynomial expressed as a linear combination of the coefficients is performed using matrix product. More precisely, if the approximating polynomial is given by the formula <\equation*> p(x)= \+\P(x) with the underlying values of to be ,x,\,x>, then the values )>, ,x> are given by <\equation*> )>>|)>>|>>|)>>>>>=(x)>|(x)>|>|(x)>>|(x)>|(x)>|>|(x)>>|>|>||>>|(x)>|(x)>|>|(x)>>>>>>>|>>|>>|>>>>>+\>|>|>>|>>>>. The t> coefficient matrix in this formula is returned by the function when called with one argument (the vector ,x,\,x)>). A call to returns a similar matrix, but with an arbitrary collection of points (not necessarily the underlying >'s). This allows us to interpolate the approximating polynomial on a fine grid. <\input| >> co = coef(yield.aov);\ <\input| >> x0 = as.numeric(as.vector(PlantDensity)); x1 = as.numeric(levels(PlantDensity)); range.x = range(x1); range.y = range(Yield); step.x = 1 <\input| >> x = seq(range.x[1], range.x[2], step.x) <\input| >> poly.coefs = poly(x1, degree = t - 1) <\input| >> pred = predict(poly.coefs, x); y.full = cbind(1, pred) %*% co; y.quadratic = cbind(1, pred) %*% c(co[1:3], 0, 0); y.cubic = cbind(1, pred) %*% c(co[1:4], 0) <\input| >> plot(range.x, range.y, type = "n", xlab = "Plant Density", ylab = "Yield", main = "The 4-th and 2-nd degree approximations") <\input| >> points(x0, Yield, col = "orange", pch = 19); points(x0, fitted(yield.aov), type = "p", col = "magenta", pch = 19); points(x, y.full, type = "l", col = "green", pch = 24); points(x, y.quadratic, type = "l", col = "blue", pch = 25); points(x, y.cubic, type = "l", col = "brown", pch = 25) <\input| >> legend("bottomright", col=c("orange","magenta","green","blue","brown"), lty = c(-1,-1,1,1,1), pch = c(19,19,-1,-1,-1), c("sample", "fitted", "4-th degree", "2-nd degree","3-rd degree")); v() <\output> |ps>||||||> <\input| >> \; > <\initial> <\collection> <\references> <\collection> > > > > > > > > > > > <\auxiliary> <\collection> <\associate|toc> |math-font-series||Data definition> |.>>>>|> |math-font-series||A plot of the data> |.>>>>|> |Setting up T|E>||||0.5fn|>|0fn|-0.1fn>>XACS||||0.5fn|>|0fn|-0.1fn>> plotting parameters |.>>>>|> > |Setting up X Windows plotting parameters |.>>>>|> > |The plot |.>>>>|> > |math-font-series||Analysis of variance> |.>>>>|> |math-font-series||Fast sums of squares for contrasts and residuals> |.>>>>|> |math-font-series||F-statistic for contrasts and the omnibus> |.>>>>|> |math-font-series||F-tests for each contrast - significance levels> |.>>>>|> |math-font-series||Check with ordinary ANOVA> |.>>>>|> |math-font-series||Plots> |.>>>>|>