Lecture 4—Wednesday, September 5, 2012

Topics

• Graphical displays of analysis of variance models
• Obtaining parameter estimates, standard errors, and confidence intervals
• Graphical display of the summary table
  • Effects graphs using the lattice package
  • Effects graphs using base graphics
• Obtaining the treatment means and their confidence intervals
  • Using the effects package
  • Using the predict function
  • Using matrix multiplication
• Interaction plots with confidence intervals
• R code

R functions and commands demonstrated

• abline is a low-level graphics command that adds vertical lines (v= ), horizontal lines (h= ), or regression lines to a plot.
• apply is used to apply functions to the rows or columns of a matrix.
• arrows draws an arrow between a pair of points on a plot.
• as.numeric converts its argument to numeric data. We used it to obtain the numeric levels of a factor.
• axis is a low-level graphics function for customizing the features of the axes of the currently displayed graph. Its first argument, 1, 2, 3, or 4, specifies the location of the axis: bottom, left, top, or right.
• bold when used in combination with expression can be used to display text in bold face
• box is a low-level graphics command that draws a box around a plot.
• coef extracts the regression coefficients from an lm object.
• confint calculates confidence intervals for the regression parameters of an lm object.
• data.frame constructs a data frame from a collection of vectors of the same length (or matrices with the same number of rows).
• diag extracts the diagonal entries from a matrix.
• dotchart is a base graphics function for producing dot plots.
• dotplot is a graphics function from the lattice package for producing dot plots.
• expression defines an R expression and is typically used to add mathematical text expressions to graphs.
• function is a keyword that begins the definition of a user-defined function.
• grid is a low-level graphics command that adds a grid of horizontal and/or vertical lines to a plot.
• legend is a low-level graphics command that adds a legend to a base graphics plot.
• lines is a low-level graphics command that draws lines between plotted points on the currently active plot.
• list concatenates objects of diverse types together as a single object.
• mtext is a low-level graphics command that can be used to add text to any one of the four margins of a graph.
• names extracts the names of the elements comprising an object.
• panel.abline (from lattice) is the panel function that draws vertical lines (v= ), horizontal lines (h= ), or regression lines.
• panel.dotplot (from lattice) is the panel function that produces a dot plot.
• panel.segments (from lattice) is the panel function that connects a pair of points with a line segment.
• panel.xyplot (from lattice) is the panel function that plots points.
• par sets global graphics parameters for base graphics.
• points is a low-level graphics command that adds individual points to the currently active plot.
• predict is used to obtain predictions (means) from a regression model for specified values of the predictors. It can also return standard errors and confidence intervals for those predictions.
• qnorm is the quantile function of a normal distribution.
• qt is the quantile function of a t-distribution.
• range calculates the minimum and maximum values of its argument and returns the result as a vector.
• rbind appends one data frame or matrix to another. The rows of the second data frame are placed after the rows of the first data frame. The column names in the two data frames must match.
• segments connects a pair of points with a line segment.
• sqrt is the square root function.
• t is the transpose function. It swaps the rows and columns of a matrix.
• vcov when applied to a regression model (lm object) returns the variance-covariance matrix of the parameter estimates.
• %*% is the matrix multiplication operator in R.

R function options

• angle= (argument to arrows function) specifies the angle of arrowhead edge; angle=90 is the appropriate setting for error bars.
• at= (argument to axis) is used to specify the locations of tick marks.
• axes= (argument to plot and other base graphics functions) can be used to prevent the x-axis and y-axis with their default settings from being drawn (by specifying axes=F).
• cex= (argument to many graphics functions) specifies the amount of character expansion to be used for plotting symbols, cex=1 is the default.
• code= (argument to arrows function) specifies the type of arrow to draw; code=3 places arrowheads on both ends.
• given.values= (argument to the effect function of the effects package) is used to specify values of model predictors that aren't part of the specified effect. Without this effect uses the mean of the remaining predictors.
• interval= (argument to predict function) specifies the type of prediction interval desired, e.g., interval="confidence".
• labels= (argument to axis and factor) specifies a vector of values to appear at the tick mark locations specified by at= in axis OR the labels for the levels defined by the factor function.
• las= (argument for many graphics text functions) controls the orientation of printed text. We chose las=2 to display tick mark labels perpendicular to the axis on which they occur.
• lend= (argument to par function) specifies the terminal ends of line segments; lend=1 yields butted ends.
• length= (argument to arrows function) specifies the size of the arrowhead.
• lineend= (argument to panel.segments function of lattice) specifies the endpoint of line segments; lineend=1 yields butted ends.
• lty= (argument to many plot functions) specifies the type of line to draw. Specifying lty=1 yields a solid line, lty=2 is a dashed line, and lty=3 is a dotted line.
• mar= (argument to par function) specifies the width of the four outer margins of plots in line units.
• mfrow= (argument to par function) specifies the arrangement of graphs in the graphics window for base graphics. For instance, mfrow=c(1,2) arranges the graphs in one row and two columns.
• newdata= (argument to predict function) specifies the values of the predictors at which to obtain predictions of a regression model.
• panel= (argument to xyplot) used to define a panel function that describes what is to be plotted in each panel.
• pch= stands for "print character" and is used to designate the plotting symbol to be used in various plotting functions: plot, points, etc.
• pt.cex= defines the size of plotting symbols in legends (overriding the value of cex).
• scales= (argument to xyplot and other lattice functions) defines characteristics of the axis—the number of tick marks, position and labels of the ticks, etc.
• se.fit= (argument to predict function) when set to TRUE causes predict to return the standard error of the predictions.
• type= (argument to most plotting functions) indicates the type of graphical display to produce. We used type='n' (where 'n' stands for 'nothing) to suppress the plotting of any data. Another option is type='b' to plot both points and line segments connecting those points.
• xlab= (argument of all plotting functions) specifies the label for the x-axis.
• xlim= (argument of all plotting functions) a vector of length 2 that specifies the minimum and maximum values to display on the x-axis.
• ylab= (argument of all plotting functions) specifies the label for the y-axis.

Additional R packages used

• effects for the effect function and its plot.eff method to generate plots of interaction effects.
• lattice for panel graphics functions.

Graphical displays of analysis of variance models

I reload the tadpoles data set and refit the analysis of variance models from last time.

We will use the second of these models, out2, to illustrate how to generate graphical displays of ANOVA models. R has three major graphics packages: base graphics, lattice, and ggplot2. The last one is a relatively new R package that is based on Leland Wilkinson's (of SYSTAT fame) "grammar of graphics". Today we focus on using base graphics and lattice to summarize ANOVA models. The plots we consider are graphs of effect estimates with error bars and interaction plots with error bars.

Obtaining parameter estimates, standard errors, and confidence intervals

ANOVA models are regression models in which all the regressors are dummy variables. This means that the regression coefficients are multiplying terms that either have value zero or one. Hence all the regression coefficients are measured on the scale and either represent means or deviations from means. Thus it makes sense to present the coefficient estimates together in the same graph. When confidence intervals are also included the magnitudes of the different effects can be compared and their relative importance assessed.

The coefficient estimates of a regression model can be extracted with the coef function.

The intercept is the mean of the reference group while the rest of the terms correspond to effects, i.e., deviations from means. To ensure that we're only comparing effects I remove the intercept from the list of estimates to plot.

To obtain the standard errors of these estimates we can extract them from the summary table of the model. They are found in column 2 of the $coefficients component of the summary object.

As was the case with the estimates I remove the intercept from the vector of standard errors.

Alternatively we can extract the standard errors from the variance-covariance matrix of the parameter estimates obtained with the vcov function. The variances lie on the diagonal which we can extract with the diag function. The standard errors are then obtained by taking their square roots with the sqrt function.

Using the estimates and their standard errors I next obtain 95% confidence intervals of the effects. Technically we should use a t-distribution for this, but the sample size used in our analysis is so large that the relative difference between using a t-distribution or a normal distribution is trivial. The degrees of freedom for the t-distribution are obtained in the $df.residual component of the model object.

To construct 95% confidence intervals we need the .025 and .975 quantiles of a t-distribution with 231 degrees of freedom. The t quantile function in R is qt.

These are barely different from the corresponding normal quantiles obtained with qnorm.

A (1 – α)% confidence interval is obtained in the usual fashion.

R provides easier ways of obtaining confidence intervals than calculating them by hand. The confint function when applied to an lm object returns confidence intervals (95% by default) for all the regression parameters in the model.

These match what we calculated by hand.

Finally I assemble everything in a data frame. We will need a column of labels for the effects. These can be obtained with the names function which extracts the labels from the vector of coefficient estimates. At the same time I make the result a factor and specify the desire3d order of the factor levels explicitly with the levels argument of factor. This prevents the data.frame function of R from putting the levels in alphabetical order.

Graphical display of the summary table

To display and compare simple point estimates such as means we can use a dot plot. A bare-bones dot plot can be obtained with the dotchart function of base graphics or the dotplot function of lattice. The syntax is slightly different for the two functions as shown below.

(a)	(b)
Fig. 1  Default dot plots as generated by (a) the dotchart function of base graphics and (b) the dotplot function of lattice

We will want to add error bars to these displays. I illustrate how to do this for lattice and base graphics separately.

Effects graphs using the lattice package

Lattice graphs are generated with a single line of code. Each lattice plotting function has a default panel function that identifies the features that will be displayed in the graph. If you want additional features beyond the default you will need to construct your own panel function. The panel function must identify everything that is to be plotted including what you would otherwise get by default. When using your own panel function it doesn't really matter what basic lattice plot function you use because the default panel function for that plot function is ignored anyway.

Panel functions can be created with the keyword panel=function(x,y) followed by a pair of curly braces { }. Special lattice panel graphing functions are then listed on separate lines within the curly braces. Each panel graphing function adds a specific feature to the graph. Typically the names of these functions begin with the word panel followed by the operation they are to perform, e.g., panel.xyplot, panel.lines, panel.points, etc., although recently lattice has been moving away from this convention.

The following lattice call uses a panel function to plot the point estimates with panel.xyplot, add a vertical line at 0 with panel.abline, and add the 95% confidence intervals with panel.segments. The xlim argument is needed to make enough room for the error bars. For this I use the minimum value of the lower 95% intervals and the maximum of the 95% confidence intervals to define the plot range. I subtract and add .02 to these values to increase this range a little bit. The x and y that appear in the declaration of the panel function correspond to the x and y-variables in the original dotplot call, in this case ests and var.labels respectively.

• panel.xyplot plots the points using filled circles, pch=16. The argument cex stands for "character expansion" and is the magnification ratio for changing the size of what's being plotted. The choice cex=1 is the default size.
• panel.abline draws a vertical line at 0, v=0, using a dashed line, lty=2, colored red, col=2.
• The syntax for panel.segments is panel.segments(x1,y1,x2,y2) where (x1, y1) and (x2, y2) are the two endpoints of the line segment to be drawn. Here the y-coordinate corresponds to the factor variable var.labels. Because the panel.segments function expects a number, not a factor, to plot, I convert the factor to a number with the as.numeric function, as.numeric(y). Fig. 2 shows the result.

Fig. 2  Dot plot with error bars (95% confidence intervals)

I redo the graph with the following improvements.

1. Lattice uses a default color for the points that turns out to be col="#0080ff" in hex notation. We can make the confidence intervals this same color by including this as an argument to the panel.segments function.
2. Although it's not obvious from Fig. 2, the default way that R draws line segments is to place a round cap at the end. The presence of the cap causes the line segment to extend beyond the coordinates that we specified. This is a problem with confidence intervals if we want to use them for hypothesis testing by determining if they cross the zero line. To override this behavior we can add the argument lineend=1 in panel.segments to obtained "butted" endpoints: line segments without caps.
3. As another cosmetic improvement we can add grid lines like those of Fig. 1. The panel.dotplot function adds grid lines by default so I could use this function first and let the rest of the function calls cover up what it plots. I color the points that panel.dotplot produces 'white' and make them small.
4. I also improve the labeling on the y-axis. I replace the hybrid names created by R that combines the factor name and the level with names that just specify the levels. I also take advantage of the mathematical typesetting abilities of R to include the multiplication symbol of arithmetic, ×, in the interaction terms. This is denoted %*% in R and requires the use of the expression function. To see other examples of mathematical typesetting in R enter demo(plotmath) at the command prompt in the console window. The R code and the mathematical expression it produces are shown in the graphics window.

To use these labels in the dotplot function we need to include the scales argument. The basic format of scales is

where x and y denote the x- and y-axes and we specify the features we want to modify on the respective axes inside the list function separated by commas. For instance, the following will cause the y-axis tick marks to be labeled using the new labels. By not specifying the x-axis I keep the default settings.

Here's the complete function call.

Fig. 3  An "improved" version of Fig. 2

So what do we learn from Fig. 3?

1. The effects graph is a graphical representation of the coefficient summary table for the model (minus the intercept term). Using Fig. 3 we can visually carry out significance tests at α = .05 by determining if the displayed confidence intervals intersect the red dashed vertical line at 0.
2. We can use Fig. 3 to assess the relative importance of the various effects in the model. This is especially useful if an appropriate reference group was chosen. In the current model the reference group is hormone treatment = corticosterone, diet = detritus, and sibship = 1. It probably would have made more sense to have made hormone treatment = normal the reference group because it serves as the control in this experiment. Putting that aside for the moment we can observe the following from Fig. 3.
   1. The hormone treatment alone produces the biggest change in the response. Both normal and Ru486 have roughly the same effect relative to the reference group corticosterone. Because normal is the control group it tells us that Ru486 has essentially no effect while corticosterone has a very large negative effect.
   2. There is no diet effect unless it is combined with Ru486. In this case switching to the shrimp diet has a weak positive effect, enough of an effect to make the Ru486 treatment significantly different from the corticosterone group.
   3. For tadpoles given a shrimp diet the increased response that is obtained by switching to Ru486 is of the same order of magnitude as the effect gained by switching from sibship 1 to sibship 2. This implies that the diet × Ru486 effect is no larger than effects attributable to ordinary background genetic variation.
3. If desired, we can also use the graph to estimate the effect attributable to any treatment combination of interest just by adding the corresponding displayed effects. For example, for a tadpole with sibship=1, diet=detritus, and hormone=Ru486 the effect is given by the single labeled estimate Ru486. For a tadpole with sibship=1, diet=shrimp, hormone=Ru486 we would need to add together the displayed effects for Ru486, Shrimp, and Ru486 × Shrimp. If we want the treatment mean instead of the effect we have to add the effect just obtained to the model intercept value which is 3.38.

Effects graphs using base graphics

With base graphics a minimal plot is created using a higher-level graphics command; we then add to that plot using lower level graphics commands. the basic strategy is to build the graph incrementally with multiple function calls. This contrasts with lattice where we create the entire graph all at once. The individual function calls in base graphics correspond exactly to the different panel functions we used in generating the lattice graph, so almost everything we did there carries over.

The lower level graphics functions of base graphics that we will use require the numeric values of the seven effect categories for plotting, so I create a new variable that numbers them 1 through 7.

When the goal is to build a modular graph, I generally use the first higher-level graphics command only to set up axis limits and labels, while generating the rest of the graph using lower-level graphics commands. Here I use the plot function to set up the y-limits and x-limits explicitly and define labels for the axes.

• In the first argument I use the range function to return the minimum and maximum x-limits (using the numerical columns of the data frame). The second argument sets the limits for the y-axis. I set these limits so that they are just beyond the plotted endpoints of 1 and 7 in order to add a little buffer around the display.

range(new.data[,2:4])
[1] -0.1380433 0.6617067
c(1-.3,7+.3)
[1] 0.7 7.3

• The result is that we have only two pairs of points to plot. By specifying type='n' , 'n' for nothing, I suppress the display of any points, so we end up only defining the x- and y-limits of the plot.
• The ylab='' argument, with nothing between the two single quotes, suppresses the printing of a label on the y-axis.
• The axes=F argument prevents the display of the x- and the y-axes. I include it because I want to add my own labels to the y-axis.

The axis functions are used to add axes to an existing plot. The first argument of axis is the desired axis. The axes are numbered 1 (bottom), 2 (left), 3 (top), and 4 (right). If you want to use default settings for the axis just enter the number and nothing more. To specify locations and labels for tick marks, add the at= and labels= argument just as was done in the scales argument of lattice graphs.

I specify las=2 to rotate the labels so that they are perpendicular to the axis. The cex.axis=.8 option reduces the size of the labels to 80% of the default size. When we look at the graph this generates we see there is a problem (Fig. 4).

Fig. 4  Graph with default margin settings

The left margin is not wide enough to accommodate the labels so they are being truncated. We need to expand the left margin beyond its default value. Margins are set with the mar argument of the par function.

The reported values are in line units and are given in the order bottom, left, top, right. I try doubling the left margin, but first I save the original settings so they can be restored later.

The box function with no arguments connects the axes and draws a box around the graph. Generating the rest of the graph is easy. I take the various panel functions from the previous section and remove the lead word "panel.". The segments function draws the confidence intervals, points adds the point estimate, and abline draws a vertical line at zero. Where an x or a y appear as the arguments in the lattice panel functions I replace them with the actual variable names. To obtain the "butted" ends on line segments requires the use of the par function again this time to set the option lend=1.

Fig. 5 Effects graph using base graphics

We can add grid lines, if desired, with the grid function. The call needed here is

The arguments nx and ny denote the number of grid lines perpendicular to the x- and y-axes respectively. If we specify nx=NA we get no grid lines. If we specify nx=NULL we get grid lines at the tick marks. Below I redo the entire graph drawing the grid lines first so they don't overwrite what has already been plotted. I then reset the margins to their default values.

Fig. 6  Effects graph with grid lines using base graphics

Obtaining the treatment means and their confidence intervals

Using the effects package

The treatment means based from an analysis of variance model can be extracted with the effect function from the effects package. In the call below I request the treatment means needed to display the fac1:fac2 interaction. In the first call the fac3 variable is set at its first level, fac32=0, using the given.values argument. In the second run I request the same treatment means at the second level of fac3. Notice that given.values expects a value for the dummy variable in the model using the name that appears in the summary table of the model, not a value for the original factor. Therefore I specify fac32=0 and fac32=1 rather than fac3=1 and fac3=2.

The component $x identifies the levels of the factors that were used to obtain the estimated means, standard errors, and confidence interval endpoints.

The components $fit, $se, $lower, and $upper contain the estimates, standard errors, and lower and upper endpoints of the 95% confidence intervals.

I assemble the output from the two effect calls in separate data frames, one for each level of fac3. In each data frame I add one more variable to indicate the level of fac3, fac3=1 or fac3=2. Finally I use the rbind function (bind by rows) to append one data frame to the other.

Using the predict function

The predict function can also be used to obtain estimates from regression models. I start by creating a data frame that contains the desired values of the predictors by using the expand.grid function to generate observations representing the levels of fac1, fac2, and fac3 in all possible combinations.

The first argument of the predict function is the regression model. The newdata argument specifies the data frame that contains the observations whose predictions we desire. I add the additional arguments se.fit=T to get standard errors and interval="confidence" to get confidence intervals. Because fac3 is currently a numeric variable in the newdata data frame, the predict function returns an error message.

I declare fac3 to be a factor and try predict again.

The $fit component of the predict object contains the estimate of the mean and the confidence interval. The $se.fit component contains the standard errors.

I assemble the results in a data frame and then rename columns 4 through 7 to match the names we used above.

Using matrix multiplication

The regression equation can be viewed as a vector dot product between a vector of values for the dummy variables and a vector of coefficient estimates. The dummy variables in the current model are defined as follows.

, , ,

Using these our final regression model out2 with two two-factor interactions can be written as follows.

where the vector c is a vector of zeros and ones and β is the parameter vector whose estimates are obtained with coef(out2). I write a function, myvec, to calculate the vector c as a function of the values of fac1, fac2, and fac3.

By plugging in values for fac1, fac2, and fac3 I obtain different values for the vector c.

To obtain treatment means I form the dot product between the vector c and the coefficient vector coef(out2). The dot product is obtained via matrix multiplication in R using the matrix multiplication operator %*%.

The matrix fac.vals2 that was created above contains all the possible combinations of the values of fac1, fac2, and fac3.

To obtain all possible c vectors I just need to evaluate the function myvec on each row of this matrix separately. The apply function of R is an easy way to do this. The apply function takes three arguments.

1. First argument: the matrix to act on
2. Second argument: the number 1 or 2. The number 1 signifies that we want to act on the rows of the matrix. The number 2 signifies that we want to act on the columns of the matrix.
3. Third argument: the function that should be applied to the rows or columns of the matrix.The function needs to be defined so that it has a single argument.

The function myvec currently takes three arguments. We need to rewrite it so that it has a single argument x with components labeled x[1], x[2], x[3]. I do this on the fly using the generic function, function(x) myvec(x[1],x[2],x[3]) when I invoke the apply function.

The results of applying myvec to the different rows of fac.vals2 appear as the columns in the output of apply. To arrange the result in the same manner as the original matrix, I transpose it using R's t function.

To obtain the treatment means for all eight treatments we form the matrix product Cβ where C is the matrix cmat just created.

Algebraically it turns out that the variance-covariance matrix of the means can be obtained with the matrix product , where Σ_β is the variance covariance matrix of the parameter estimates obtained with vcov(out2). The second appearance of C in the product is the matrix transpose of the first. To obtain the standard errors of the means we extract the diagonal of the final matrix and take its square root.

The 95% confidence intervals for the mean can now be obtained in the usual fashion in which we use the residual degrees of freedom from the model for the degrees of freedom of the t-distribution.

I assemble the results in a data frame.

Interaction plots with confidence intervals

I illustrate plotting the fac1 × fac2 interaction separately for different levels of fac3. Each graph will consist of the mean profiles of fac1 separately by the levels of fac2 for a single value of fac3. This yields two graphs, one for each value of fac3. The goal is to produce the following two graphs, here generated using the effects package, but with confidence intervals displayed for each of the means.

(a)	(b)
Fig. 7  Interaction plots of the fac1:fac2 interaction as produced by the effects package

Because fac1 is to appear on the x-axis I need the numerical values of the three factor levels. I obtain them by using the as.numeric function on the factor variable.

As I did with the effects graph, I start by setting up the graph but plotting nothing. I use the minimum of the low95 values and the maximum of the up95 values to establish the y-limits. The plotted x-values are numbered 1, 2, 3, so I set the x-limits to begin just a little before 1 and to end just a little after 3. I then use the axis function to add labels to the tick marks on the x-axis. I use the levels function to extract the labels from the factor variable fac1.

Next I subset the rows of the data frame so that we first plot only the means corresponding to fac3 = 1.

I want to draw the mean profile for each diet (fac2) group separately. So I subset the data again by the first diet, fac2 = 'D'.

To prevent the error bars of the two mean profiles from lying on top of each other I shift one profile a little bit to the left and the other profile a little bit to the right. The variable myjitter contains the amount of the shift for the first profile.

The arrows function is used to draw the confidence intervals. Its basic syntax is the same as the segments function from before but it supports the additional arguments angle, code, and length.

• The angle argument controls the shape of the arrowhead; angle=90 yields flat ends.
• The code argument determines where the arrowheads should appear; code=3 puts arrowheads on both ends.
• The length argument controls the length of the arrowhead. This will typically need adjustment but for these data length=.05 works well.

I follow this up with the lines function to draw the mean profile and the points function to plot the estimates at the middle of the confidence intervals. They each use the same syntax, x-variable followed by the y-variable. Just like the confidence intervals each has to shifted the same amount.

I next repeat these same lines of code for the second diet, fac2 = 'S', with the following changes.

• I change the subsetting level so I select observations for which fac2 = 'S'.
• I use a positive value of myjitter so that the points and error bars are moved slightly to the right.
• I change the color of the profile to red, col=2, and change the symbol type to open circles, pch=1. First I plot the points using filled white circles, col='white' and pch=16, to cover up the drawn lines. Then I add red open circles on top.

Finally I add a small header on top of the graph with the mtext (for margin text) function to indicate the level of fac3 that was plotted, fac3 = 1.

Lastly I add a legend to identify which diet corresponds to which mean profile.

• The first argument of legend is the location of the legend in the plot. The nine allowable key words that identify the position are "bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center". It is also possible to specify the location of the legend by explicitly giving the x- and y-coordinates of the top left corner of the legend.
• The second argument of legend is the vector of text labels that are to appear in the legend.
• I specify the remaining arguments by name. They include the colors, col, and symbols, pch, used in the plot. These are listed in an order that corresponds to the order of the text labels.
• The cex argument controls the size of text and symbols in the legend. It is possible to make the text and symbols different sizes in the legend. The argument pt.cex can be used to control the symbol size in which case cex controls the size of everything else.
• I add a title to the legend with the title argument. I make the title bold face by using a mathematical expression and the bold function.
• The argument bty='n' suppresses the drawing of a box around the legend.

Fig. 8  Mean profile plot for sibship 1

To plot the mean profiles for sibship 2 I make only three changes to the code listed above.

1. I subset the original data frame so sibship 2 (fac3=2) is selected instead of sibship 1.
2. I change the header on top of the graph to indicate that sibship 2 was plotted.
3. I drop the legend because it already appears in the first graph.

Notice that in the initial plot function I use all of the data to set the y-limits, not just the data for sibship 2. This ensures that the scale is the same in both graphs.

Fig. 9  Mean profile plot for sibship 2

Finally I display both plots in the same graphics window side by side. For this I use the mfrow argument of par to divide the graphics window into one row and two columns.

Fig. 10  Mean profile plots for the tadpole experiment

Notice that although we have explicitly plotted only the fac1 × fac2 interaction, both interactions are clearly on display. The fac1 × fac2 interaction is obvious because in each panel the red and black profiles are not parallel, but we also can see evidence of the fac2 × fac3 interaction. If fac3 was a purely additive effect then the right panel would be identical to the left panel except for the fact that the profiles would be shifted up or down (the main effect of fac3). Instead what we see is that the two diet profiles (fac2) are further apart in the right panel (sibship 2) than in the left panel (sibship 1). Thus sibship (fac3) is modifying the effect of diet. It has modified it sufficiently that now in the left panel the diet effect is significantly different from zero at all three hormonal treatments (the confidence intervals fail to overlap).

R Code

A compact collection of all the R code displayed in this document appears here.

Course Home Page