R Tutorial Series

R Tutorial Series: Centering Variables and Generating Z-Scores with the Scale() Function

noreply@blogger.com (John Quick) — Thu, 01 Mar 2012 14:00:00 +0000

Centering variables and creating z-scores are two common data analysis activities. While they are relatively simple to calculate by hand, R makes these operations extremely easy thanks to the scale() function.

Tutorial Files

Before we begin, you may want to download the dataset (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory.

The Scale() Function

The scale() function makes use of the following arguments.

x: a numeric object
center: if TRUE, the objects' column means are subtracted from the values in those columns (ignoring NAs); if FALSE, centering is not performed
scale: if TRUE, the centered column values are divided by the column's standard deviation (when center is also TRUE; otherwise, the root mean square is used); if FALSE, scaling is not performed

Centering Variables

Normally, to center a variable, you would subtract the mean of all data points from each individual data point. With scale(), this can be accomplished in one simple call.

> #center variable A using the scale() function

> scale(A, center = TRUE, scale = FALSE)

You can verify these results by making the calculation by hand, as demonstrated in the following screenshot.

Centering a variable with the scale() function and by hand

Generating Z-Scores

Normally, to create z-scores (standardized scores) from a variable, you would subtract the mean of all data points from each individual data point, then divide those points by the standard deviation of all points. Again, this can be accomplished in one call using scale().

> #generate z-scores for variable A using the scale() function

> scale(A, center = TRUE, scale = TRUE)

Again, the following screenshot demonstrates equivalence between the function results and hand calculation.

Generating z-scores from a variable by hand and using the scale() function

Complete Scale() Example

To see a complete example of how scale() can be used to center variables and generate z-scores in R, please download the scale() example (.txt) file.

References

The official scale function manual page is available from: http://stat.ethz.ch/R-manual/R-patched/library/base/html/scale.html

R Tutorial Series: Citing R with EndNote

noreply@blogger.com (John Quick) — Fri, 06 Jan 2012 15:55:00 +0000

Unfortunately, due to the vexing complexities of academic style guides and the limitations of associated software packages, citing a non-standard name, such as Cher, Prince, or R Development Core Team can be problematic. Thankfully, I have discovered a simple trick in Word and EndNote that allows for the accurate automatic formatting of R citations. Note that this method was developed using Word 2011 and EndNote X4 for Mac. I am unaware of the differences between operating systems and software versions, but it is anticipated that this method will work for almost anyone.

The Intuitive, But Nonworking Way

If you were going to create your R record in EndNote, you would probably enter something like what is pictured below. In the name field, it makes sense to just type in R Development Core Team.

However, this is where things take an untimely turn. EndNote will try to interpret that peculiar name as a series of first, last, and middle names, which leads to inaccurate citations.

The Unintiuitive, But Working Way

This is where we basically need to trick EndNote into interpreting our R citation the proper way. All we have to do is add a comma after R Development Core Team in the name field.

This tells EndNote that R Core Development Team is a complete last name of an author that has no first name. Hence, EndNote uses what it has (a last name with no first name) in generating its citations.

Note: The official citation for R can be found by issuing the citation() command in the R console.

R Tutorial Series: Exploratory Factor Analysis

noreply@blogger.com (John Quick) — Mon, 24 Oct 2011 14:00:00 +0000

Exploratory factor analysis (EFA) is a common technique in the social sciences for explaining the variance between several measured variables as a smaller set of latent variables. EFA is often used to consolidate survey data by revealing the groupings (factors) that underly individual questions. This will be the context for demonstration in this tutorial.

Tutorial Files

Before we begin, you may want to download the dataset (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 300 responses on 6 items from a survey of college students' favorite subject matter. The items range in value from 1 to 5, which represent a scale from Strongly Dislike to Strongly Like. Our 6 items asked students to rate their liking of different college subject matter areas, including biology (BIO), geology (GEO), chemistry (CHEM), algebra (ALG), calculus (CALC), and statistics (STAT). This is where our tutorial ends, because all students rated all of these content areas as Strongly Dislike, thereby rendering insufficient variance for conducting EFA (just kidding).

Beginning Steps

To begin, we need to read our datasets into R and store their contents in variables.

> #read the dataset into R variable using the read.csv(file) function

> data <- read.csv("dataset_EFA.csv")

First 10 rows of the dataset

Psych Package

Next, we need to install and load the psych package, which I prefer to use when conducting EFA. In this tutorial, we will make use of the package's fa() function.

> #install the package

> install.packages("psych")

> #load the package

> library(psych)

Number of Factors

For this tutorial, we will assume that the appropriate number of factors has already been determined to be 2, such as through eigenvalues, scree tests, and a priori considerations. Most often, you will want to test solutions above and below the determined amount to ensure the optimal number of factors was selected.

Factor Solution

To derive the factor solution, we will use the fa() function from the psych package, which receives the following primary arguments.

r: the correlation matrix
nfactors: number of factors to be extracted (default = 1)
rotate: one of several matrix rotation methods, such as "varimax" or "oblimin"
fm: one of several factoring methods, such as "pa" (principal axis) or "ml" (maximum likelihood)

Note that several rotation and factoring methods are available when conducting EFA. Rotation methods can be described as orthogonal, which do not allow the resulting factors to be correlated, and oblique, which do allow the resulting factors to be correlated. Factoring methods can be described as common, which are used when the goal is to better describe data, and component, which are used when the goal is to reduce the amount of data. The fa() function is used for common factoring. For component analysis, see princomp(). The best methods will vary by circumstance and it is therefore recommended that you seek professional council in determining the optimal parameters for your future EFAs.
In this tutorial, we will use oblique rotation (rotate = "oblimin"), which recognizes that there is likely to be some correlation between students' latent subject matter preference factors in the real world. We will use principal axis factoring (fm = "pa"), because we are most interested in identifying the underlying constructs in the data.

> #calculate the correlation matrix

> corMat <- cor(data)

> #display the correlation matrix

> corMat

The correlation matrix

> #use fa() to conduct an oblique principal-axis exploratory factor analysis

> #save the solution to an R variable

> solution <- fa(r = corMat, nfactors = 2, rotate = "oblimin", fm = "pa")

> #display the solution output

> solution

Complete solution output

By looking at our factor loadings, we can begin to assess our factor solution. We can see that BIO, GEO, and CHEM all have high factor loadings around 0.8 on the first factor (PA1). Therefore, we might call this factor Science and consider it representative of a student's interest in science subject matter. Similarly, ALG, CALC, and STAT load highly on the second factor (PA2), which we might call Math. Note that STAT has a much lower loading on PA2 than ALG or CALC and that it has a slight loading on factor PA1. This suggests that statistics is less related to the concept of Math than algebra and calculus. Just below the loadings table, we can see that each factor accounted for around 30% of the variance in responses, leading to a factor solution that accounted for 66% of the total variance in students' subject matter preference. Lastly, notice that our factors are correlated at 0.21 and recall that our choice of oblique rotation allowed for the recognition of this relationship.
Of course, there are many other considerations to be made in developing and assessing an EFA that will not be presented here. The intent with this tutorial was simply to demonstrate the basic execution of EFA in R. For a detailed and digestible overview of EFA, I recommend the Factor Analysis chapter of Multivariate Data Analysis by Hair, Black, Babin, and Anderson.

Complete EFA Example

To see a complete example of how EFA data can be organized using the psych package in R, please download the EFA example (.txt) file. For the code used in this tutorial, download the EFA Example (.R) file.

References

Revelle, W. (2011). psych: Procedures for Personality and Psychological Research. http://personality-project.org/r/psych.manual.pdf

The R Programming Wikibook

noreply@blogger.com (John Quick) — Fri, 08 Jul 2011 14:00:00 +0000

The R Programming wikibook is an open source community project that "aims to create a cross-disciplinary practical guide to the R programming language." It was launched in June 2011 and is seeking content and contributors. The full call for the R Programming wikibook can be found on Tal Galili's blog. The R Programming wikibook itself is available at http://en.wikibooks.org/wiki/R_Programming.

I am writing to raise awareness for the R Programming wikibook and to formally offer content from the R Tutorial Series by John M. Quick for use under a Creative Commons Attribution-ShareAlike 3.0 Unported License. This means that articles from the R Tutorial Series may be included in and modified for the R Programming wikibook, so long as proper attribution is given and the resulting content is made available under an equivalent license. A complete list of contributing blogs can be found on the R Programming wikibook's Sources page. I hope that the R Programming wikibook will thrive and grow to support a large community of R users, including readers of the R Tutorial Series.

R Tutorial Series: 2011 ANOVA Article Data

noreply@blogger.com (John Quick) — Mon, 28 Mar 2011 14:00:00 +0000

Having wrapped up a recent flurry of R ANOVA articles (and exhausted my knowledge of the subject), I decided to take a look at the R Tutorial Series' Google Analytics data from the past few months.
Since I posted the Two-Way Omnibus ANOVA article on January 17, we have had about 150 visits per day and over 19,000 total page views. The original introduction to R posts are still the most popular ones here, although a few from the regression and ANOVA series are also represented in the most viewed.
I also wanted to share a funny observation about the patterns of visits to the R Tutorial Series. The following graph portrays our daily viewership over the past few months.

The valleys on the chart correspond to the weekends (evidently when no one wants to read about statistical computing). The initial peaks are Mondays, which happen to be when I most often make new posts. Typically, a slightly higher peak comes on Tuesday, followed by a gradual decline back into the weekend valley. Thus, our visits to the R Tutorial Series end up creating a nice little wave pattern throughout the year. The good news with these numbers is that people are reading the R Tutorial Series and (hopefully) learning to use R and apply it to their daily work, which is the blog's ultimate purpose. I also appreciate all of the comments, questions, and tips that have been posted by readers. Your feedback really helps to improve the tutorials.

Upcoming Plans

As mentioned, I have concluded my planned coverage of ANOVA in R. Thus, we have reached a sort of break period, much like the one that followed last year's spurt of regression tutorials. I do plan to keep writing R tutorials, but for the time being, they may arrive in less predictable intervals and cover a wider variety of content. As always, I welcome contact regarding guest posts, especially on statistical and R content that I have not covered or am not yet familiar with.

R Tutorial Series: Applying the Reshape Package to Organize ANOVA Data

noreply@blogger.com (John Quick) — Mon, 14 Mar 2011 14:00:00 +0000

As demonstrated in the preceding ANOVA tutorials, data organization is central to conducting ANOVA in R. In standard ANOVA, we used the tapply() function to generate a table for a single summary function. In repeated measures ANOVA, we used separate datasets for our omnibus ANOVA and follow-up comparisons. This tutorial will demonstrate how the reshape package can be used to simplify the ANOVA data organization process in R.

Tutorial Files

Before we begin, you may want to download the between group and repeated measures datasets (.csv) used in this tutorial. Be sure to right-click and save the files to your R working directory. The between groups dataset contains a hypothetical sample of 30 cases separated into three groups (a, b, and c). The repeated measures dataset contains a hypothetical sample of 10 cases across three measurements (a, b, and c). In both cases, the values are represented on a scale that ranges from 1 to 5.

Beginning Steps

To begin, we need to read our datasets into R and store their contents in variables.

> #read the datasets into R variables using the read.csv(file) function

> dataBetween <- read.csv("dataset_ANOVA_reshape_1.csv")

> dataRepeated <- read.csv("dataset_ANOVA_reshape_2.csv")

Reshape Package

Next, we need to install and load the reshape package. In this tutorial, we will make use of the package's cast() and melt() functions.

> #install the package

> install.packages("reshape")

> #load the package

> library(reshape)

Using cast() to Derive ANOVA Descriptives

The cast() function can be used to easily derive summary statistics for a between groups ANOVA dataset. The cast() function receives the following primary arguments.

data: the dataset
formula: in our case, a one-sided formula indicating the grouping variable
fun.aggregate: a function or vector of functions for deriving summary statistics, such as mean, var, or sd

> #display the raw between groups data

> dataBetween

The raw between groups data

> #cast the between groups data using cast(data, formula, fun.aggregate) to get the group means

> cast(dataBetween, formula = ~group, fun.aggregate = mean)

The casted data with means

Note that the fun.aggregate argument can also receive a vector of summary statistics functions. This will yield all of the requested descriptives via a single cast() function.

> #cast the between groups data using cast(data, formula, fun.aggregate) to get the group means, variances, and standard deviations

> cast(dataBetween, formula = ~group, fun.aggregate = c(mean, var, sd))

The casted data with descriptives

Using melt() to Prepare Repeated Measures Data for Pairwise Comparisons

The melt() function can be used to morph a repeated measures ANOVA dataset prior to conducting pairwise comparisons. The melt() function receives the following primary arguments.

data: the dataset
id.vars: the id variable or a vector of values that can be used as ids
measure.vars: a vector containing the variables to be melted
variable_name: the name of the column containing the melted variables

> #display the repeated measures data

> dataRepeated

The raw repeated measures data

> #melt the repeated measures data using melt(data, id.vars, measure.vars, variable_name) to organize it for pairwise comparisons

> melt(dataRepeated, id.vars = "case", measure.vars = c("valueA", "valueB", "valueC"), variable_name = "abcValues")

The melted repeated measures data

Note that the data are now prepared to be used in the pairwise.t.test() function. See the One-Way ANOVA with Pairwise Comparisons tutorial for details on using the pairwise.t.test() function.

Complete ANOVA Reshape Example

To see a complete example of how ANOVA data can be organized using the reshape package in R, please download the ANOVA reshape example (.txt) file.

R Tutorial Series: ANOVA Pairwise Comparison Methods

noreply@blogger.com (John Quick) — Mon, 07 Mar 2011 14:00:00 +0000

When we have a statistically significant effect in ANOVA and an independent variable of more than two levels, we typically want to make follow-up comparisons. There are numerous methods for making pairwise comparisons and this tutorial will demonstrate how to execute several different techniques in R.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 30 participants who are divided into three stress reduction treatment groups (mental, physical, and medical). The values are represented on a scale that ranges from 1 to 5. This dataset can be conceptualized as a comparison between three stress treatment programs, one using mental methods, one using physical training, and one using medication. The values represent how effective the treatment programs were at reducing participant's stress levels, with higher numbers indicating higher effectiveness.

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the dataset into an R variable using the read.csv(file) function

> dataPairwiseComparisons <- read.csv("dataset_ANOVA_OneWayComparisons.csv")

> #display the data

> dataPairwiseComparisons

The first ten rows of our dataset

Omnibus ANOVA

For the purposes of this tutorial, we will assume that the omnibus ANOVA has already been conducted and that the main effect for treatment was statistically significant. For details on this process, see the One-Way ANOVA with Pairwise Comparisons tutorial, which uses the same dataset.

Means

Let's also look at the means of our treatment groups. Here, we will use the tapply() function, along with the following arguments, to generate a table of means.

X: the data
INDEX: a list() of factor variables
FUN: the function to be applied

> #use tapply(X, INDEX, FUN) to generate a table displaying each treatment group mean

> tapply(X = dataPairwiseComparisons$StressReduction, INDEX = list(dataPairwiseComparisons$Treatment), FUN = mean)

The treatment group means

Pairwise Comparisons

We will cover five major techniques for controlling Type I error when making pairwise comparisons. These methods are no adjustment, Bonferroni's adjustment, Holm's adjustment, Fisher's LSD, and Tukey's HSD. All of these techniques will be demonstrated on our sample dataset, although the decision as to which to use in a given situation is left up to the reader.

pairwise.t.test()

Our first three methods will make use of the pairwise.t.test() function, which has the following major arguments.

x: the dependent variable
g: the independent variable
p.adj: the p-value adjustment method used to control for the family-wise Type I error rate across the comparisons; one of "none", "bonferroni", "holm", "hochberg", "hommel", "BH", or "BY"

No Adjustment

Using p.adj = "none" in the pairwise.t.test() function makes no correction for the Type I error rate across the pairwise tests. This technique can be useful for employing methods that are not already built into R functions, such as the Shaffer/Modified Shaffer, which use different alpha level divisors based on the number of levels composing the independent variable. The console results will contain no adjustment, but the researcher can manually consider the statistical significance of the p-values under his or her desired alpha level.

> #use pairwise.t.test(x, g, p.adj) to test the pairwise comparisons between the treatment group means

> #no adjustment

> pairwise.t.test(dataPairwiseComparisons$StressReduction, dataPairwiseComparisons$Treatment, p.adj = "none")

Pairwise comparisons of treatment group means with no adjustment

With no adjustment, the mental-medical and physical-medical comparisons are statistically significant, whereas the mental-physical comparison is not. This suggests that both the mental and physical treatments are superior to the medical treatment, but that there is insufficient statistical support to distinguish between the mental and physical treatments.

Bonferroni Adjustment

The Bonferroni adjustment simply divides the Type I error rate (.05) by the number of tests (in this case, three). Hence, this method is often considered overly conservative. The Bonferroni adjustment can be made using p.adj = "bonferroni" in the pairwise.t.test() function.

> #Bonferroni adjustment

> pairwise.t.test(dataPairwiseComparisons$StressReduction, dataPairwiseComparisons$Treatment, p.adj = "bonferroni")

Pairwise comparisons of treatment group means using Bonferroni adjustment

Using the Bonferroni adjustment, only the mental-medical comparison is statistically significant. This suggests that the mental treatment is superior to the medical treatment, but that there is insufficient statistical support to distinguish between the mental and physical treatments and the physical and medical treatments. Notice that these results are more conservative than with no adjustment.

Holm Adjustment

The Holm adjustment sequentially compares the lowest p-value with a Type I error rate that is reduced for each consecutive test. In our case, this means that our first p-value is tested at the .05/3 level (.017), second at the .05/2 level (.025), and third at the .05/1 level (.05). This method is generally considered superior to the Bonferroni adjustment and can be employed using p.adj = "holm" in the pairwise.t.test() function.

> #Holm adjustment

> pairwise.t.test(dataPairwiseComparisons$StressReduction, dataPairwiseComparisons$Treatment, p.adj = "holm")

Pairwise comparisons of treatment group means using Holm adjustment

Using the Holm procedure, our results are practically (but not mathematically) identical to using no adjustment.

LSD Method

The Fisher Least Significant Difference (LSD) method essentially does not correct for the Type I error rate for multiple comparisons and is generally not recommended relative to other options. However, should the need arise to employ this method, one should seek out the LSD.test() function in the agricolae package, which has the following major arguments.

y: the dependent variable
trt: the independent variable
DFerror: the degrees of freedom error
MSerror: the mean squared error

Note that the DFerror and MSerror can be found in the omnibus ANOVA table.

> #load the agricolae package (install first, if necessary)

> library(agricolae)

#LSD method

#use LSD.test(y, trt, DFerror, MSerror) to test the pairwise comparisons between the treatment group means

> LSD.test(dataPairwiseComparisons$StressReduction, dataPairwiseComparisons$Treatment, 30.5, 1.13)

Pairwise comparisons of treatment group means using LSD method

Using the LSD method, our results are practically (but not mathematically) identical to using no adjustment or the Holm procedure.

HSD Method

The Tukey Honest Significant Difference (HSD) method controls for the Type I error rate across multiple comparisons and is generally considered an acceptable technique. This method can be executed using the TukeyHSD(x) function, where x is a linear model object created using the aov(formula, data) function. Note that in this application, the aov(formula, data) function is identical to the lm(formula, data) that we are already familiar with from linear regression.

> #HSD method

> #use TukeyHSD(x), in tandem with aov(formula, data), to test the pairwise comparisons between the treatment group means

TukeyHSD(aov(StressReduction ~ Treatment, dataPairwiseComparisons))

Pairwise comparisons of treatment group means using HSD method

Using the HSD method, our results are practically (but not mathematically) identical to using the Bonferroni, Holm, or LSD methods.

Complete Pairwise Comparisons Example

To see a complete example of how various pairwise comparison techniques can be applied in R, please download the ANOVA pairwise comparisons example (.txt) file.

R Tutorial Series: Two-Way ANOVA with Unequal Sample Sizes

noreply@blogger.com (John Quick) — Mon, 28 Feb 2011 14:00:00 +0000

When the sample sizes within the levels of our independent variables are not equal, we have to handle our ANOVA differently than in the typical two-way case. This tutorial will demonstrate how to conduct a two-way ANOVA in R when the sample sizes within each level of the independent variables are not the same.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 30 students who were exposed to one of two learning environments (offline or online) and one of two methods of instruction (classroom or tutor), then tested on a math assessment. Possible math scores range from 0 to 100 and indicate how well each student performed on the math assessment. Each student participated in either an offline or online learning environment and received either classroom instruction (i.e. one to many) or instruction from a personal tutor (i.e. one to one).

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the dataset into an R variable using the read.csv(file) function

> dataTwoWayUnequalSample <- read.csv("dataset_ANOVA_TwoWayUnequalSample.csv")

> #display the data

> dataTwoWayUnequalSample

The first ten rows of our dataset

Unequal Sample Sizes

In our study, 16 students participated in the online environment, whereas only 14 participated in the offline environment. Further, 20 students received classroom instruction, whereas only 10 received personal tutor instruction. As such, we should take action to compensate for the unequal sample sizes in order to retain the validity of our analysis. Generally, this comes down to examining the correlation between the factors and the causes of the unequal sample sizes en route to choosing whether to use weighted or unweighted means - a decision which can drastically impact the results of an ANOVA. This tutorial will demonstrate how to conduct ANOVA using both weighted and unweighted means. Thus, the ultimate decision as to the use of weighted or unweighted means is left up to each individual and his or her specific circumstances.

Weighted Means

First, let's suppose that we decided to go with weighted means, which take into account the correlation between our factors that results from having treatment groups with different sample sizes. A weighted mean is calculated by simply adding up all of the values and dividing by the total number of values. Consequently, we can easily derive the weighted means for each treatment group using our subset(data, condition) and mean(data) functions.

> #use subset(data, condition) to create subsets for each treatment group

> #offline subset

> offlineData <- subset(dataTwoWayUnequalSample, dataTwoWayUnequalSample$environment == "offline")

> #online subset

> onlineData <- subset(dataTwoWayUnequalSample, dataTwoWayUnequalSample$environment == "online")

> #classroom subset

> classroomData <- subset(dataTwoWayUnequalSample, dataTwoWayUnequalSample$instruction == "classroom")

> #tutor subset

> tutorData <- subset(dataTwoWayUnequalSample, dataTwoWayUnequalSample$instruction == "tutor")

> #use mean(data) to calculate the weighted means for each treatment group

> #offline weighted mean

> mean(offlineData$math)

> #online weighted mean

> mean(onlineData$math)

> #classroom weighted mean

> mean(classroomData$math)

> #tutor weighted mean

> mean(tutorData$math)

The weighted means for the environment and instruction conditions

ANOVA using Type I Sums of Squares

When applying weighted means, it is suggested that we use Type I sums of squares (SS) in our ANOVA. Type I happens to be the default SS used in our standard anova(object) function, which will be used to execute our analysis. Note that in the case of two-way ANOVA, the ordering of our independent variables matters when using weighted means. Therefore, we must run our ANOVA two times, once with each independent variable taking the lead. However, the interaction effect is not affected by the ordering of the independent variables.

> #use anova(object) to execute the Type I SS ANOVAs

> #environment ANOVA

> anova(lm(math ~ environment * instruction, dataTwoWayUnequalSample))

> #instruction ANOVA

> anova(lm(math ~ instruction * environment, dataTwoWayUnequalSample))

The Type I SS ANOVA results. Note the differences in main effects based on the ordering of the independent variables.

These results indicate statistically insignificant main effects for both the environment and instruction variables, as well as the interaction between them.

Unweighted Means

Now let's turn to using unweighted means, which essentially ignore the correlation between the independent variables that arise from unequal sample sizes. An unweighted mean is calculated by taking the average of the individual group means. Thus, we can derive our unweighted means by summing the means of each level of our independent variables and dividing by the total number of levels. For instance, to find the unweighted mean for environment, we will add the means for our offline and online groups, then divide by two.

> #use mean(data) and subset(data, condition) to calculate the unweighted means for each treatment group

> #offline unweighted mean = (classroom offline mean + tutor offline mean) / 2

(mean(subset(offlineData$math, offlineData$instruction == "classroom")) + mean(subset(offlineData$math, offlineData$instruction == "tutor"))) / 2

> #online unweighted mean = (classroom online mean + tutor online mean) / 2

> (mean(subset(onlineData$math, onlineData$instruction == "classroom")) + mean(subset(onlineData$math, onlineData$instruction == "tutor"))) / 2

> #classroom unweighted mean = (offline classroom mean + online classroom mean) / 2

> (mean(subset(classroomData$math, classroomData$environment == "offline")) + mean(subset(classroomData$math, classroomData$environment == "online"))) / 2

> #tutor unweighted mean = (offline tutor mean + online tutor mean) / 2

> (mean(subset(tutorData$math, tutorData$environment == "offline")) + mean(subset(tutorData$math, tutorData$environment == "online"))) / 2

The unweighted means for the environment and instruction conditions

ANOVA using Type III Sums of Squares

When applying unweighted means, it is suggested that we use Type III sums of squares (SS) in our ANOVA. Type III SS can be set using the type argument in the Anova(mod, type) function, which is a member of the car package.

> #load the car package (install first, if necessary)

> library(car)

> #use the Anova(mod, type) function to conduct the Type III SS ANOVA

> Anova(lm(math ~ environment * instruction, dataTwoWayUnequalSample), type = "3")

The Type III SS ANOVA results.

Once again, our ANOVA results indicate statistically insignificant main effects for both the environment and instruction variables, as well as the interaction between them. However, it is worth noting that both the means and p-values are different when using unweighted means and Type III SS compared to weighted means and Type I SS. In certain cases, this difference can be quite pronounced and lead to entirely different outcomes between the two methods. Hence, choosing the appropriate means and SS for a given analysis is a matter that should be approached with conscious consideration.

Pairwise Comparisons

Note that since our independent variables contain only two levels, there is no need to conduct follow-up comparisons. However, should you reach this point with a statistically significant independent variable of more than three levels, you could conduct pairwise comparisons in the same manner as demonstrated in the Two-Way ANOVA with Comparisons tutorial.

Complete Two-Way ANOVA with Unequal Sample Sizes Example

To see a complete example of how two-way ANOVA with unequal sample sizes can be conducted in R, please download the two-way ANOVA with unequal sample sizes example (.txt) file.

R Tutorial Series: Two-Way Repeated Measures ANOVA

noreply@blogger.com (John Quick) — Mon, 21 Feb 2011 14:00:00 +0000

Repeated measures data require a different analysis procedure than our typical two-way ANOVA and subsequently follow a different R process. This tutorial will demonstrate how to conduct two-way repeated measures ANOVA in R using the Anova() function from the car package.
Note that the two-way repeated measures ANOVA process can be very complex to organize and execute in R. Although it has been distilled into just a few small steps in this guide, it is recommended that you fully and precisely complete the example before experimenting with your own data. As you will see, organization of the raw data is critical to successfully conducting a two-way repeated measures ANOVA using the demonstrated technique.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) and sample idata frame (.csv) used in this tutorial. Be sure to right-click and save the files to your R working directory. This dataset contains a hypothetical sample of 30 participants whose interest in school and interest in work was measured at three different ages (10, 15, and 20). The interest values are represented on a scale that ranges from 1 to 5 and indicate how interested each participant was in a given topic at each given age.

Data Setup

Notice that our data are arranged differently for a repeated measures ANOVA. In a typical two-way ANOVA, we would place all of the values of our independent variable in a single column and identify their respective levels with a second column, as demonstrated in this sample two-way dataset. In a two-way repeated measures ANOVA, we instead combine each independent variable with its time interval, thus yielding columns for each pairing. Hence, rather than having one vertical column for school interest and one for work interest, with a second column for age, we have six separate columns for interest, three for school interest and three for work interest at each age level. The following graphic is intended to help demonstrate this organization method.

Treat time as if it were an independent variable. Then combine each independent variable with each level of time and arrange the columns horizontally.

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the dataset into an R variable using the read.csv(file) function

> dataTwoWayRepeatedMeasures <- read.csv("dataset_ANOVA_TwoWayRepeatedMeasures.csv")

> #display the data

> #notice the atypical column arrangement for repeated measures data

> dataTwoWayRepeatedMeasures

The first ten rows of our dataset

idata Frame

Another item that we need to import for this analysis is our idata frame. This object will be used in our Anova() function to define the structure of our analysis.

> #read the idata frame into an R variable

> idataTwoWayRepeatedMeasures <- read.csv("idata_ANOVA_TwoWayRepeatedMeasures.csv")

> #display the idata frame

> #notice the text values and correspondence between our idata rows and the columns in our dataset

> idataTwoWayRepeatedMeasures

The idata frame

Note that it is critical that your idata frame take the demonstrated form for this technique to work. I experimented with several alternative, perhaps more intuitive, layouts without success. It is particularly important to notice that both columns of the idata frame contain text values (not numerical ones - hence the repeated prefixing of Age to the values in every row of the Age column). Additionally, if you read the rows of the idata frame horizontally, you will see that they correspond precisely to the columns of our dataset. The following graphic is intended to help demonstrate this organization method.

Use only text values in your idata frame. Ensure that the rows of your idata frame correspond to the columns in your dataset.

Linear Model

Prior to executing our analysis, we must follow two steps to formulate our linear model to be used in the Anova() function.

Step 1: Bind the Columns

> #use cbind() to bind the columns of the original dataset

> interestBind <- cbind(dataTwoWayRepeatedMeasures$schoolAge10, dataTwoWayRepeatedMeasures$schoolAge15, dataTwoWayRepeatedMeasures$schoolAge20, dataTwoWayRepeatedMeasures$workAge10, dataTwoWayRepeatedMeasures$workAge15, dataTwoWayRepeatedMeasures$workAge20)

Step 2: Define the Model

> #use lm() to generate a linear model using the bound columns from step 1

> interestModel <- lm(interestBind ~ 1)

Anova(mod, idata, idesign) Function

Typically, researchers will choose one of several techniques for analyzing repeated measures data, such as an epsilon-correction method, like Huynh-Feldt or Greenhouse-Geisser, or a multivariate method, like Wilks' Lambda or Hotelling's Trace. Conveniently, having already prepared our data, we can employ a single Anova(mod, idata, idesign) function from the car package to yield all of the relevant repeated measures results. This allows us simplicity in that only a single function is required, regardless of the technique that we wish to employ. Thus, witnessing our outcomes becomes as simple as locating the desired method in the cleanly printed results.
Our Anova(mod, idata, idesign) function will be composed of three arguments. First, mod will contain our linear model. Second, idata will contain our data frame. Third, idesign will contain a multiplication of the row headings from our idata frame (in other words, our independent variables), preceded by a tilde (~). Thus, our final function takes on the following form.

> #load the car package (install first, if necessary)

library(car)

> #compose the Anova(mod, idata, idesign) function

> analysis <- Anova(interestModel, idata = idataTwoWayRepeatedMeasures, idesign = ~Interest * Age)

Results Summary

Finally, we can use the summary(object) function to visualize the results of our repeated measures ANOVA.

> #use summary(object) to visualize the results of the repeated measures ANOVA

> summary(analysis)

Relevant segment of repeated measures ANOVA results

Supposing that we are interested in the Wilks' Lambda method, we can see that there is a statistically significant interaction effect between interest in school and interest in work across the age groups (p < .001). This suggests that we should further examine our data at the level of simple main effects. For more information investigating on simple main effects, see the Two-Way ANOVA with Interactions and Simple Main Effects tutorial. Of course, in this case of repeated measures ANOVA, another way to break the data down would be to run two one-way repeated measures ANOVAs, one for each of the independent variables. In either instance, pairwise comparisons can be conducted to determine the significance of the differences between the levels of any significant effects.

Complete Two-Way Repeated Measures ANOVA Example

To see a complete example of how two-way repeated measures ANOVA can be conducted in R, please download the two-way repeated measures ANOVA example (.txt) file.

References

Moore, Colleen. (n.d.). 610 R9 -- Two-way Repeated-measures Anova. Retrieved January 21, 2011 from http://psych.wisc.edu/moore/Rpdf/610-R9_Within2way.pdf

Book Review: R Graphs Cookbook

noreply@blogger.com (John Quick) — Thu, 17 Feb 2011 14:00:00 +0000

Book Information

Mittal, H. (2011). R graphs cookbook. Birmingham, UK: Packt Publishing Ltd.

Audience

The book's stated audience is anyone who is familiar with the basics of R, as well as expert users who are looking for a graphical reference. However, it is my opinion that the book is better suited for advanced users who are already somewhat familiar with R graphics and are very comfortable with programming in R.

Content

To begin, the first chapter of R Graphs Cookbook rapidly introduces all of the major graphic types covered in the book. Next, in Chapter two, readers are acquainted with various arguments and modification functions that are used throughout the book to customize and enhance visuals. Subsequently, individual chapters focus on specific topics in R graphics, such as:

scatterplots,
line and time series charts,
bar, dot, and pie charts,
histograms,
box plots,
heat and contour maps,
geographical maps, and
exporting and annotating graphics.

Analysis

I will start with some general impressions, before moving into chapter by chapter analyses.

First, I feel that the book needs both more and larger screenshots. Often times, recipes are without any visuals and most of the time only one is present, whereas one per major graphical modification is expected. Furthermore, the screenshots are too small. These are critical items to neglect in a book that explicitly deals with visuals. Fortunately, full-size, full-color images are provided with the downloadable code for the book.

Second, I feel that the topics presented in the book are glanced over with far too little explanation. This is the main reason that I feel it is not suited for those who are not already well versed in R programming. Moreover, R Graphs Cookbook frequently refers the reader to help documentation or to other books on R, which can be frustrating. I personally feel that a book should be largely self-contained, at least when discussing topics within its scope.

Third, I believe that the book could be better organized for use as a fast reference guide and that it generally could be better structured to present information. For example, rather than tables are a clearer way to present head to head comparisons between objects, and lists are better for describing several function arguments.

On the other hand, I do like the book's code formatting, which displays one argument per line. While this could confuse novice users into thinking that each argument is a separate line of executable code, most readers should find this a welcomed organization style for often lengthy graphics functions. I also enjoyed how the see also sections at the end of each recipe let me know whether more recipes would build on a given topic.

Continuing, chapter one felt like a whirlwind of information that charged forward with a lack of purpose, organization, and explanation. Chapter two was much better, offering several nice recipes that were fast and easy to digest, with just enough information provided.

Chapter three takes an in-depth look at scatterplots and provides a number of useful recipes, such as how to group data, label points, generate error bars, and create graphical correlation matrices. Similarly, chapter four provides a solid collection of recipes for time series and line charts.

In contrast, chapters five through seven cover a disappointingly sparse amount of material related to their respective topics. Unfortunately, they do not stretch far beyond what is covered in the two graphics-focused chapters of Statistical Analysis with R, which is a guide for newcomers and early beginners. From an advanced reference like R Graphs Cookbook, I expected broader coverage. For instance, very few external packages are presented in this book, with the author choosing to focus on built-in graphics functions almost exclusively. An introduction to external options, such as ggplot2, would be warmly welcomed.

Chapters eight and nine relate a few of the lesser covered topics in R, including heat, contour, and geographical maps. These chapters will likely be informative and valuable for readers interested in these graphical applications.

Lastly, chapter ten deals with the presentation and exportation of graphics. While I wish a deeper exploration was made, there are some useful tips in is chapter. Namely, the use of the expression() function to annotate graphics is well covered.

Brief Summary

Title: R Graphs Cookbook
Author: Hrishi Mittal
Where To Find: Packt Publishing
Audience: those who are comfortable programming in R, able to mix, match, apply, and extend recipes for their own purposes, and looking to learn more about R's built-in graphical capabilities.
Content: a loosely associated collection of recipes for applying R's built-in graphics functions to create the most common types of charts, graphs, plots, and maps.
Analysis: although it could have better visuals, structure, and coverage, it is likely that almost any reader will be able to take away valuable techniques from this book
Arbitrary Rating: 6/10
Recommendation: take a look at the table of contents and count the number of recipes that would both be useful to you and that you do not already know how to accomplish to get an idea of how much you will take away from this book; also read the free sample chapter
Disclaimer: I received a review copy of this book

R Tutorial Series: One-Way Repeated Measures ANOVA

noreply@blogger.com (John Quick) — Mon, 14 Feb 2011 14:00:00 +0000

Repeated measures data require a different analysis procedure than our typical one-way ANOVA and subsequently follow a different R process. This tutorial will demonstrate how to conduct one-way repeated measures ANOVA in R using the Anova(mod, idata, idesign) function from the car package.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 30 participants whose interest in voting was measured at three different ages (10, 15, and 20). The interest values are represented on a scale that ranges from 1 to 5 and indicate how interested each participant was in voting at each given age.

Data Setup

Notice that our data are arranged differently for a repeated measures ANOVA. In a typical one-way ANOVA, we would place all of the values of our independent variable in a single column and identify their respective levels with a second column, as demonstrated in this sample one-way dataset. In a repeated measures ANOVA, we instead treat each level of our independent variable as if it were a variable, thus placing them side by side as columns. Hence, rather than having one vertical column for voting interest, with a second column for age, we have three separate columns for voting interest, one for each age level.

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the dataset into an R variable using the read.csv(file) function

> dataOneWayRepeatedMeasures <- read.csv("dataset_ANOVA_OneWayRepeatedMeasures.csv")

> #display the data

> #notice the atypical column arrangement for repeated measures data

> dataOneWayRepeatedMeasures

The first ten rows of our dataset

Preparing the Repeated Measures Factor

Prior to executing our analysis, we must follow a small series of steps in order to prepare our repeated measures factor.

Step 1: Define the Levels

> #use c() to create a vector containing the number of levels within the repeated measures factor

> #create a vector numbering the levels for our three voting interest age groups

> ageLevels <- c(1, 2, 3)

Step 2: Define the Factor

> #use as.factor() to create a factor using the level vector from step 1

> #convert the age levels into a factor

> ageFactor <- as.factor(ageLevels)

Step 3: Define the Frame

> #use data.frame() to create a data frame using the factor from step 2

> #convert the age factor into a data frame

> ageFrame <- data.frame(ageFactor)

Step 4: Bind the Columns

> #use cbind() to bind the levels of the factor from the original dataset

> #bind the age columns

> ageBind <- cbind(dataOneWayRepeatedMeasures$Interest10, dataOneWayRepeatedMeasures$Interest15, dataOneWayRepeatedMeasures$Interest20)

Step 5: Define the Model

> #use lm() to generate a linear model using the bound factor levels from step 4

> #generate a linear model using the bound age levels

> ageModel <- lm(ageBind ~ 1)

Employing the Anova(mod, idata, idesign) Function

Typically, researchers will choose one of several techniques for analyzing repeated measures data, such as an epsilon-correction method, like Huynh-Feldt or Greenhouse-Geisser, or a multivariate method, like Wilks' Lambda or Hotelling's Trace. Conveniently, having already prepared our data, we can employ a single Anova(mod, idata, idesign) function from the car package to yield all of the relevant repeated measures results. This allows us simplicity in that only a single function is required, regardless of the technique that we wish to employ. Thus, witnessing our outcomes becomes as simple as locating the desired method in the cleanly printed results.
Our Anova(mod, idata, idesign) function will be composed of three arguments. First, mod will contain our linear model from Step 5 in the preceding section. Second, idata will contain our data frame from Step 3. Third, idesign will contain our factor from Step 2, preceded by a tilde (~). Thus, our final function takes on the following form.

> #load the car package (install first, if necessary)

library(car)

> #compose the Anova(mod, idata, idesign) function

> analysis <- Anova(ageModel, idata = ageFrame, idesign = ~ageFactor)

Visualizing the Results

Finally, we can use the summary(object) function to visualize the results of our repeated measures ANOVA.

> #use summary(object) to visualize the results of the repeated measures ANOVA

> summary(analysis)

Relevant segment of repeated measures ANOVA results

Supposing that we are interested in the Wilks' Lambda method, we can see that the differences in the means for voting interest at ages 10, 15, and 20 are statistically significant (p < .001).

Pairwise Comparisons

Note that we could conduct follow-up comparisons on our age factor to determine which age level means are significantly different from one another. For this purpose, it is recommended that the data be rearranged into the standard ANOVA format that we have used throughout our other tutorials. Subsequently, we could conduct pairwise comparisons in the same manner as demonstrated in the One-Way ANOVA with Comparisons tutorial.

Complete One-Way Repeated Measures ANOVA Example

To see a complete example of how one-way repeated measures ANOVA can be conducted in R, please download the one-way repeated measures ANOVA example (.txt) file.

R Tutorial Series: Two-Way ANOVA with Interactions and Simple Main Effects

noreply@blogger.com (John Quick) — Mon, 07 Feb 2011 14:00:00 +0000

When an interaction is present in a two-way ANOVA, we typically choose to ignore the main effects and elect to investigate the simple main effects when making pairwise comparisons. This tutorial will demonstrate how to conduct pairwise comparisons when an interaction is present in a two-way ANOVA.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 60 participants who are divided into three stress reduction treatment groups (mental, physical, and medical) and two gender groups (male and female). The stress reduction values are represented on a scale that ranges from 1 to 5. This dataset can be conceptualized as a comparison between three stress treatment programs, one using mental methods, one using physical training, and one using medication across genders. The values represent how effective the treatment programs were at reducing participant's stress levels, with higher numbers indicating higher effectiveness.

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the dataset into an R variable using the read.csv(file) function

> dataTwoWayInteraction <- read.csv("dataset_ANOVA_TwoWayInteraction.csv")

> #display the data

> dataTwoWayInteraction

The first ten rows of our dataset.

Omnibus Test

Let's run a general omnibus test to assess the main effects and interactions present in the dataset.

> #use anova(object) to test the omnibus hypothesis

> #Are main effects or interaction effects present in the independent variables?

> anova(lm(StressReduction ~ Treatment * Gender, dataTwoWayInteraction))

The omnibus ANOVA test

Divide the Data

The significant omnibus interaction suggests that we should ignore the main effects and instead investigate the simple main effects for our independent variables. To do so, we need to divide our dataset along each level of our treatment variable. We can create subsets of our dataset using the subset(data, condition) function, where data is the original dataset and condition contains the parameters defining the subset.

> #use subset(data, condition) to divide the original dataset

> #medical subset

> dataMedical <- subset(dataTwoWayInteraction, Treatment == "medical")

> #mental subset

> dataMental <- subset(dataTwoWayInteraction, Treatment == "mental")

> #physical subset

> dataPhysical <- subset(dataTwoWayInteraction, Treatment == "physical")

Group ANOVAs

With datasets representing each of our treatment groups, we can now run an ANOVA for each that investigates the impact of gender. You may notice that this is effectively running three one-way ANOVAs with gender being the independent variable. Therefore, we should control for Type I error by dividing our typical .05 alpha level by three (.017).

> #run ANOVA on the treatment subsets to investigate the impacts of gender within each

> #medical

> anova(lm(StressReduction ~ Gender, dataMedical))

> #mental

> anova(lm(StressReduction ~ Gender, dataMental))

> #physical

> anova(lm(StressReduction ~ Gender, dataPhysical))

The gender within treatment group ANOVA tests

At an alpha level of .017, the gender effect within the mental (p = .014) and physical (p < .001) groups was statistically significant. In the mental condition, the means are 3 for males and 4 for females. In the physical condition, the means are 4 for males and 2 for females. These results suggest that the mental treatment is more effective in reducing stress for females than males, while the physical treatment is more effective for males than females. Further, there is insufficient statistical support for a gender difference in the medical treatment.

Pairwise Comparisons

Note that since our gender variable contains only two levels, there is no need to conduct follow-up comparisons. However, should you reach this point with an independent variable of more than three levels, you could conduct pairwise comparisons in the same manner as demonstrated in the Two-Way ANOVA with Comparisons tutorial. In this case, remember to carry through your reduced Type I error rate from the preceding ANOVA tests.

Complete Two-Way ANOVA with Interactions Example

To see a complete example of how two-way ANOVA simple main effects can be explored in R, please download the two-way ANOVA interaction example (.txt) file.

Statistical Analysis with R Book Reviews

noreply@blogger.com (John Quick) — Wed, 02 Feb 2011 14:00:00 +0000

Reviews of my Statistical Analysis with R book have started to emerge online and I am writing today to share them with potential readers and recommenders.

Reviews

The following is a list of online reviews for Statistical Analysis with R. If you have written a review of the book and would like it to be featured in this post, please contact me.

In summarizing the reviews, a few points are very clear about Statistical Analysis with R.

It is for beginners: The book was written for people who have little to no experience with R, statistical software, and programming. It makes no assumptions of prior experience along these lines and starts right from the beginning. If you are new to R and want to learn how to apply it to your work, then this book is for you. If you are already an intermediate or experienced user, perhaps you might recommend it to people you know who are just becoming familiar with R.
It is a learning tool, not a reference: The book is structured with the intent that it is experienced as a holistic learning experience. The chapters build on one another and progressively delve deeper into R. It is not a dictionary-style reference that one might pull out, flip to an entry, and get a brief answer on a single item. Again, this has implications for the audience. Beginners are more likely to enjoy this approach, whereas experienced users may be interested in more of a reference-style book.
It has a story: Woven into the book's learning structure is a storyline based on the Three Kingdoms period of ancient Chinese history. For many, this will be a motivating and engaging way to learn. For others, the story may not inspire the same level of interest. If you would like to get a taste of the story, and the book in general, it is recommended that you read the free sample chapter.

New Release: R Graph Cookbook

Packt Publishing recently released a second book on R, the R Graph Cookbook by Hrishi Mittal. This reference-style guide covers an array of R graphical applications and is geared towards users who are already familiar with the basics of R. A sample chapter is available. I will be reviewing this book in the near future and posting my thoughts here on the R Tutorial Series blog. [Update: read my review of R Graph Cookbook]

R Tutorial Series: Two-Way ANOVA with Pairwise Comparisons

noreply@blogger.com (John Quick) — Mon, 31 Jan 2011 14:00:00 +0000

By extending our one-way ANOVA procedure, we can test the pairwise comparisons between the levels of several independent variables. This tutorial will demonstrate how to conduct pairwise comparisons in a two-way ANOVA.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 27 participants who are divided into three stress reduction treatment groups (mental, physical, and medical) and three age groups (young, mid, and old). The stress reduction values are represented on a scale that ranges from 0 to 10. This dataset can be conceptualized as a comparison between three stress treatment programs, one using mental methods, one using physical training, and one using medication, across three age groups. The stress reduction values represent how effective the treatment programs were at reducing participant's stress levels, with higher numbers indicating higher effectiveness. Note that the numbers in this dataset are not very realistic and are simply used to make this example possible.

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the dataset into an R variable using the read.csv(file) function

> dataTwoWayComparisons <- read.csv("dataset_ANOVA_TwoWayComparisons.csv")

> #display the data

> dataTwoWayComparisons

The first ten rows of our dataset.

Omnibus Test

Let's run a general omnibus test to assess the main effects and interactions present in the dataset.

> #use anova(object) to test the omnibus hypothesis

> #Are main effects or interaction effects present in the independent variables?

> anova(lm(StressReduction ~ Treatment * Age, dataTwoWayComparisons))

The omnibus ANOVA test

Pairwise Comparisons

Since the omnibus test was significant for both variables and no interaction effect was present, we can proceed to testing the main effect pairwise comparisons. To accomplish this, we will apply our pairwise.t.test() function to each of our independent variables. For more details on the pairwise.t.test() function, see the One-Way ANOVA with Pairwise Comparisons tutorial.

> #use pairwise.t.test(x, g, p.adj) to test the pairwise comparisons between the treatment group means

> #What significant differences are present amongst the treatment means?

> pairwise.t.test(dataTwoWayComparisons$StressReduction, dataTwoWayComparisons$Treatment, p.adj = "none")

> #use pairwise.t.test(x, g, p.adj) to test the pairwise comparisons between the age group means

> #What significant differences are present amongst the age group means?

> pairwise.t.test(dataTwoWayComparisons$StressReduction, dataTwoWayComparisons$Age, p.adj = "none")

Pairwise comparisons of treatment group means

Pairwise comparisons of age group means

Note that the desired p-adjustment method will vary by researcher, study, etc. Here, we will assume an alpha level of .05 for all tests, effectively making no adjustment for the family-wise Type I error rate.
These results indicate that there are are no statistically significant pairwise differences between the treatment groups and that all of the comparisons between age groups are statistically significant. The age group means are 8 for young, 5 for mid, and 2 for old. Consequently, we are inclined to conclude that, regardless of treatment, young patients are going to be most responsive, followed by middle aged patients, followed by older ones. However, there is insufficient support to differentiate between the effectiveness of the treatment methods themselves.

Complete Two-Way ANOVA with Pairwise Comparisons Example

To see a complete example of how two-way ANOVA pairwise comparisons can be conducted in R, please download the two-way ANOVA comparisons example (.txt) file.

R Tutorial Series: One-Way ANOVA with Pairwise Comparisons

noreply@blogger.com (John Quick) — Mon, 24 Jan 2011 14:00:00 +0000

When we have more than two groups in a one-way ANOVA, we typically want to statistically assess the differences between each group. Whereas a one-way omnibus ANOVA assesses whether a significant difference exists at all amongst the groups, pairwise comparisons can be used to determine which group differences are statistically significant.

Tutorial Files

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the dataset into an R variable using the read.csv(file) function

> dataOneWayComparisons <- read.csv("dataset_ANOVA_OneWayComparisons.csv")

> #display the data

> dataOneWayComparisons

The first ten rows of our dataset.

Omnibus Test

One way to begin an ANOVA is to run a general omnibus test. The advantage to starting here is that if the omnibus test comes up insignificant, you can stop your analysis and deem all pairwise comparisons insignificant. If the omnibus test is significant, you should continue with pairwise comparisons.

> #use anova(object) to test the omnibus hypothesis

> #Is there a significant difference amongst the treatment means?

> anova(lm(StressReduction ~ Treatment, dataOneWayComparisons))

The omnibus ANOVA test

Pairwise Comparisons

Since the omnibus test was significant, we are safe to continue with our pairwise comparisons. To make pairwise comparisons between the treatment groups, we will use the pairwise.t.test() function, which has the following major arguments.

x: the dependent variable
g: the independent variable
p.adj: the p-value adjustment method used to control for the family-wise Type I error rate across the comparisons; one of "none", "bonferroni", "holm", "hochberg", "hommel", "BH", or "BY"

The pairwise.t.test() function can be implemented as follows.

> #use pairwise.t.test(x, g, p.adj) to test the pairwise comparisons between the treatment group means

> #What significant differences are present amongst the treatment means?

> pairwise.t.test(dataOneWayComparisons$StressReduction, dataOneWayComparisons$Treatment, p.adj = "none")

Pairwise comparisons of treatment group means

Note that the desired p-adjustment method will vary by researcher, study, etc. Here, we will assume an alpha level of .05 for all tests, effectively making no adjustment for the family-wise Type I error rate.
These results indicate that there is a statistically significant difference between the mental and medical (p = .004) and physical and medical (p = 0.045), but not the mental and physical (p = 0.302) treatments. The treatment means are 3.5 for mental, 3 for physical, and 2 for medical. Subsequently, we are inclined to conclude based on this study that the mental and physical treatments lead to greater stress reduction than the medical method, but that there is insufficient statistical support to determine that either the mental or physical treatment method is superior.

Complete One-Way ANOVA with Pairwise Comparisons Example

To see a complete example of how one-way ANOVA pairwise comparisons can be conducted in R, please download the one-way ANOVA comparisons example (.txt) file.

References

R Documentation (n.d.). Adjust P-values for Multiple Comparisons. Retrieved January 16, 2011 from http://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html

R Tutorial Series: Two-Way Omnibus ANOVA

noreply@blogger.com (John Quick) — Mon, 17 Jan 2011 14:00:00 +0000

As with the one-way case, testing the omnibus hypothesis via two-way ANOVA is simple process in R. This tutorial will explore how R can be used to perform a two-way ANOVA to test the difference between two (or more) group means.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 60 participants who are divided by gender (male and female) and treatment group (control and treatment). The values represent a scale that ranges from 1 to 5. For instance, this dataset could be conceptualized as a comparison between two professional training programs, where the control group participated the company's longstanding program and the treatment group participated in an experimental program. The values could represent the attitudes of employees towards the training programs on a scale from 1 (poor) to 5 (excellent).

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the two-way ANOVA dataset into an R variable using the read.csv(file) function

> dataTwoWay <- read.csv("dataset_ANOVA_TwoWay.csv")

> #display the data

> dataTwoWay

The first ten rows of our two-way ANOVA dataset.

Two-Way ANOVA

Now that our data are ready, we can conduct a two-way omnibus ANOVA test using the anova(object) function. Note that the only step necessary to add a second independent variable into our ANOVA model is to incorporate it into our lm(model, dataset) function using the * operator. Whereas our one-way model was lm(Values ~ Group), our two-way model becomes lm(Values ~ Group * Gender). As you can see from the results below, adding a second independent variable in this manner also gives us information about the interaction between our variables.

> #use anova(object) to test the omnibus hypothesis in two-way ANOVA

> #Are the differences between the group means for treatment and gender statistically significant?

> #Is there a statistically significant interaction between treatment and gender?

> anova(lm(Values ~ Group * Gender, dataTwoWay))

Our two-way ANOVA table.

The output of our ANOVA test indicates that the difference between our treatment group means is statistically significant (p < .001) and that the difference between genders is not (p = .585). However, in light of the statistically significant interaction between treatment group and gender (p = .032), we would generally elect to forgo these main effects. Subsequently, a series of follow-up procedures could be carried out to examine the simple main effects for treatment group and gender.

Two-Way Multiple Group ANOVA

Conducting a two-way omnibus ANOVA with multiple groups is identical to the demonstrated two-group test. The only difference is that the values in your dataset would be associated with more than two groups. Subsequently, the omnibus hypothesis would test for mean differences across all of the groups. The anova(object) function and its contained lm(formula, data) function would remain the same.

Complete Two-Way Omnibus ANOVA Example

To see a complete example of how a two-way omnibus ANOVA can be conducted in R, please download the two-way ANOVA example (.txt) file.

R Beginner's Guide Book Update: Statistical Analysis with R Released

noreply@blogger.com (John Quick) — Mon, 08 Nov 2010 14:00:00 +0000

In the final days of October, my beginner's guide to R was released. The book's official title is Statistical Analysis with R and it can be found on the Packt Publishing website.
The primary focus of Statistical Analysis with R is helping new users become accustomed to R and empowering them to apply R to suit their own needs. No prior experience with R, statistical software packages, or programming is necessary to learn from this book. It is written for a broad audience and should be well received by businesspeople, IT professionals, researchers, and students alike. Statistical Analysis with R takes readers on a journey from their first installation and launch of R, to analyzing and assessing data, to communicating and visualizing results. This guide is an excellent way to rapidly become an experienced R user and learn the skills that you need to apply R to your work.

Samples

A sample chapter from Statistical Analysis with R is available from the Packt website. This chapter, the book's eight, introduces the graphical capabilities of R, such as generating, customizing, and exporting various plots, charts, and graphs. You can download the sample chapter and its accompanying R files for free. If you like this chapter and are interested in learning more about R's graphical capabilities, you should know that chapter 9 demonstrates in depth how you can build and customize your own R visualizations.
The publisher has also posted a few brief samples from the book, which can be accessed via the following links. These samples are taken from chapters 7 and 8 of Statistical Analysis with R. Respectively, they cover the common process behind all R analyses and introduce the graphical capabilities of R.

Feedback

If you decide to read Statistical Analysis with R, please feel free to provide me with your feedback. It would be great to know what you learned from the book, how future guides could be improved, and your overall experience with R.

R Tutorial Series: One-Way Omnibus ANOVA

noreply@blogger.com (John Quick) — Mon, 11 Oct 2010 21:01:00 +0000

Testing the omnibus hypothesis via one-way ANOVA is simple process in R. This tutorial will explore how R can be used to perform a one-way ANOVA to test the difference between two (or more) group means.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains a hypothetical sample of 60 participants, who are divided into two groups (control and treatment) of 30. The values represent a scale that ranges from 1 to 5. For instance, this dataset could be conceptualized as a comparison between two professional training programs, where the control group participated the company's longstanding program and the treatment group participated in an experimental program. The values could represent the attitudes of employees towards the training programs on a scale from 1 (poor) to 5 (excellent).

Beginning Steps

To begin, we need to read our dataset into R and store its contents in a variable.

> #read the one-way ANOVA dataset into an R variable using the read.csv(file) function

> dataOneWay <- read.csv("dataset_ANOVA_OneWay.csv")

> #display the data

> dataOneWay

The first ten rows of our one-way ANOVA dataset.

One-Way ANOVA

Now that our data are ready, we can conduct a one-way omnibus ANOVA test using the anova(object) function.

> #use anova(object) to test the omnibus hypothesis in one-way ANOVA

> #is the difference between the group means statistically significant?

> anova(lm(Values ~ Group, dataOneWay))

Our one-way ANOVA table.

The output of our ANOVA test indicates that the difference between our group means is statistically significant (p < .001). Conceptually, this suggests that employee attitudes towards the experimental training program were significantly higher than their attitudes towards the preexisting program.

Note that the object argument in our anova(object) function contained a linear model generated by the lm(formula, data) function. This is the same type of model that is used when conducting linear regression in R. A more detailed explanation of the lm(formula, data) function and examples of its use are available in my Simple Linear Regression article.

One-Way Multiple Group ANOVA

Conducting a one-way omnibus ANOVA with multiple groups is identical to the demonstrated two-group test. The only difference is that the values in your dataset would be associated with more than two groups. Subsequently, the omnibus hypothesis would test for mean differences across all of the groups. The anova(object) function and its contained lm(formula, data) function would remain the same.

Complete One-Way Omnibus ANOVA Example

To see a complete example of how a one-way omnibus ANOVA can be conducted in R, please download the one-way ANOVA example (.txt) file.

R Beginner's Guide Book Update 10/1/2010

noreply@blogger.com (John Quick) — Fri, 01 Oct 2010 14:00:00 +0000

Update: Statistical Analysis with R is now available!

I recently submitted the final drafts of all chapters of my R Beginner's Guide book, which is to be published through Packt. The official publishing timeline is set to December 2010, although the book may release ahead of schedule if all continues to go well. Below is an updated list of the major topics covered in the R Beginner's Guide.

Over the course of this book, you will acquire the knowledge and skills necessary to:

Conduct organized data analyses in R
Communicate data analyses conducted in R
Generate, customize, and export detailed charts, plots, and graphs
Build your own custom data visualizations
Program in the R language
Create your own custom functions
Extend the functionality of R via external packages
Manage the R workspace and console
Import external data into R
Manipulate data using variables
Execute a wide array of multi-argument and variable-argument functions
Develop and employ predictive regression models
Assess the practical and statistical significance of predictions
Understand R, its benefits, and how to use it to maximize the impact of your data analyses

R Tutorial Series: Labeling Data Points on a Plot

noreply@blogger.com (John Quick) — Sun, 19 Sep 2010 18:05:00 +0000

There are times that labeling a plot's data points can be very useful, such as when conveying information in certain visuals or looking for patterns in our data. Fortunately, labeling the individual data points on a plot is a relatively simple process in R. In this tutorial, we will use the Calibrate package's textxy function to label the points on a scatterplot.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains information used to estimate undergraduate enrollment at the University of New Mexico (Office of Institutional Research, 1990). Note that this tutorial assumes that this data has already been read into R and saved into a variable named enrollmentData.

Plot

To begin, we need to create a scatterplot using the plot(x,y) function. With our example data, we will plot the year on the x axis and the unemployment rate on the y axis.

> #generate a plot using the plot(x,y) function

> #plot year on the x axis and unemployment rate on the y axis

> plot(enrollmentData$YEAR, enrollmentData$UNEM)

For a more detailed description of plotting data in R, see the article on scatterplots.

Textxy

Within the calibrate package, the textxy() function can be used to label a plot's data points. The textxy() function accepts the following arugments ("Label points in a plot," n.d.).

x: the x values of the plot's points
y: the y values of the plot's points
labs: the labels to be associated with the plot's points

cx: used to resize the label font
dcol: used to set the label color; defaults to black
m: sets the origin of the plot; defaults to (0,0)

Here, we will use textxy() to add labels for the enrollment at the University of New Mexico to each of our plot's data points.

> #if necessary, install the calibrate package

> #install.packages("calibrate")

> #load the calibrate package

> library(calibrate)

> #use the textxy() function to add labels to the preexisting plot's points

> #add labels for the total enrollment

> textxy(enrollmentData$YEAR, enrollmentData$UNEM, enrollmentData$ROLL)

In this case, adding labels to our data points helps us to better assess the relationships in our dataset.

Complete Data Point Labeling Example

To see a complete example of how a plot's data points can be labeled in R, please download the Data Point Labeling (.txt) file.

References

Label points in a plot. (n.d.). Retrieved September 19, 2010 from http://rss.acs.unt.edu/Rdoc/library/calibrate/html/textxy.html
Office of Institutional Research (1990). Enrollment Forecast [Data File]. Retrieved November 22, 2009 from http://lib.stat.cmu.edu/DASL/Datafiles/enrolldat.html

R Beginner's Guide Book Update 7/19/2010

noreply@blogger.com (John Quick) — Mon, 19 Jul 2010 17:46:00 +0000

Update: Statistical Analysis with R is now available!

I am excited to announce that I have submitted the entire first draft of my R Beginner's Guide book, which is to be published through Packt. The tenth and final chapter was submitted a full month ahead of schedule. The printed book could become available in as little as three to four months.

Below is a list of the major topics covered in the R Beginner's Guide.

Understanding what R is, its benefits, and why to use it
Downloading, installing, and running R
Dissecting the anatomy of R
Programming in R
Handling external data
Using variables
Managing the R workspace and console
Using multi-argument and variable-argument functions
Creating predictive data models
Assessing practical vs. statistical significance
Regression modeling
Creating custom functions
Assessing the viability of predictions
Organizing and communicating data analyses
Generating, customizing, and exporting graphics
Building custom visualizations
Extending R via packages
Taking advantage of electronic learning resources

R Beginner's Guide Book Update 4/28/2010

noreply@blogger.com (John Quick) — Thu, 29 Apr 2010 00:44:00 +0000

Update: Statistical Analysis with R is now available!

I am writing to update you on the progress of my R Beginner's Guide book, which is to be published through Packt. I have really gotten to work over the past couple months and have recently completed the first draft of the first half of the book. Right now, I am operating a few weeks ahead of our planned schedule, which calls for the first draft of all ten chapters by mid-August.

To give you an idea of its content, the book focuses on most of the topics covered in this blog as well as many more, such as data visualization, custom functions, and online resources. The topics are covered in great depth and numerous opportunities for practice and exploration are offered. The book's theme centers around the Three Kingdoms period of ancient China. The reader takes on the role of the lead strategist for the Shu kingdom at a pivotal point in history. Throughout the book, the reader uses R to help devise a course of action for the Shu forces.

I will continue to make steady progress on this book over the summer months. I am also excited to be able to share more R tutorials and knowledge in the near future through this project.

R Tutorial Series: R Beginner's Guide and R Bloggers Updates

noreply@blogger.com (John Quick) — Mon, 15 Mar 2010 14:00:00 +0000

1/1/2011 Update: Tal Galili wrote an article that revisits the first year of R-Bloggers and this post was listed as one of the top 14. Therefore, I decided to make a small update to each section. I start by describing the initial series of tutorials that I wrote. A few more have been added since and even more planned in the upcoming year. As always, an up to date listing of my articles can be found on the R Tutorial Series blog. New posts will also continue to be offered through the R Bloggers network.

Since October 2009, I have written 13 articles [many more now, of course] for the R Tutorial Series blog. The first two introduce new users to R. The remaining 11 cover a wide range of topics related to multiple regression and correlation. This collection of tutorials represents my most recent training in statistics. Thus, for the time being, I will not be contributing new articles as frequently as I have over the past few months. However, I will undoubtedly encounter future projects that require new statistical methods and partake in more statistics courses, both of which will provide additional tutorial material. Below is a categorized list of the articles currently offered in the R Tutorial Series.

Introduction to R

Descriptive Statistics

Summary and Descriptive Statistics

Data Visualization

Scatterplots

Correlation

Zero-Order Correlations

Regression

I also have two additional R-related items to update you on. The first is the R Bloggers website and the second is my R Beginner's Guide.

1/1/2011 Update: I originally reported that 50 blogs composed the R Bloggers network. Now that number has risen to over 140. I hope that R Bloggers continues to thrive and contribute to the R community.

R Tutorial Series on R Bloggers

R Bloggers (http://www.r-bloggers.com) is a website that aggregates over 50 different blogs that focus on R. It is an excellent resource for keeping up to date on the many uses of R and for learning about the wide range of work being conducted in R. I recommend using R Bloggers for these purposes. The R Tutorial Series was invited to participate in the R Bloggers collection and is now available to R Bloggers' readers.

R Beginner's Guide

11/1/2010 Update: Statistical Analysis with R is now available!

Lastly, I want to let you know that I am working on a beginner's guide for R. It is primarily focused towards introducing R to information technology, business, and data analyst professionals. The book will be offered through PACKT Publishing (http://www.packtpub.com) and should be available within the next year. If you have enjoyed the R Tutorial Series, then you may be interested in looking for the guide once it is completed. In the meantime, keep reading the R Tutorial Series and R Bloggers and I will keep you updated on the book's major milestones.

R Tutorial Series: Basic Polynomial Regression

noreply@blogger.com (John Quick) — Mon, 08 Feb 2010 14:00:00 +0000

Often times, a scatterplot reveals a pattern that seems not so linear. Polynomial regression can be used to explore a predictor at different levels of curvilinearity. This tutorial will demonstrate how polynomial regression can be used in a hierarchical fashion to best represent a dataset in R.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. Note that all code samples in this tutorial assume that this data has already been read into an R variable and has been attached. This dataset contains hypothetical student data that uses practice exam scores to predict final exam scores.

Scatterplot

The preceding scatterplot demonstrates that these data may not be linear. Notably, no one scored lower than 50 on the practice exam and at approximately the 85 and above practice mark, final exam scores taper off. These suggest that the data is curvilinear. Furthermore, since exam scores range between 0 to 100, it is not possible to observe nor appropriate to predict that an individual with a 150 practice score would have a certain final exam score.

Creating The Higher Order Variables

A two step process, identical to the one used to create interaction variables, can be followed to create higher order variables in R. First, the variables must be centered to mitigate multicollinearity. Second, the predictor must be multiplied by itself a certain number of times to create each higher order variable. In this tutorial, we will explore the a linear, quadratic, and cubic model. Therefore, the predictor will need to be squared to create the quadratic model and cubed to create the cubic model.

Step 1: Centering

To center a variable, simply subtract its mean from each data point and save the result into a new R variable, as demonstrated below.

> #center the independent variable

> FinalC <- Final - mean(Final)

> #center the predictor

> PracticeC <- Practice - mean(Practice)

Step 2: Multiplication

Once the input variable has been centered, the higher order terms can be created. Since a higher order variable is formed by the product of a predictor with itself, we can simply multiply our centered term from step one and save the result into a new R variable, as demonstrated below.

> #create the quadratic variable

> PracticeC2 <- PracticeC * PracticeC

> #create the cubic variable

> PracticeC3 <- PracticeC * PracticeC * PracticeC

Creating The Models

Now we have all of the pieces necessary to assemble our linear and curvilinear models.

> #create the models using lm(FORMULA, DATAVAR)

> #linear model

> linearModel <- lm(FinalC ~ PracticeC, datavar)

> #quadratic model

> quadraticModel <- lm(FinalC ~ PracticeC + PracticeC2, datavar)

> #cubic model

> cubicModel <- lm(FinalC ~ PracticeC + PracticeC2 + PracticeC3, datavar)

Evaluating The Models

As is the case in other forms of regression, it can be helpful to summarize and compare our potential models using the summary(MODEL) and anova(MODEL1, MODEL2,… MODELi) functions.

> #display summary information about the models
> summary(linearModel)
> summary(quadraticModel)
> summary(cubicModel)
#compare the models using ANOVA
anova(linearModel, quadraticModel, cubicModel)

The model summaries and ANOVA comparison chart are displayed below.

At this point we can compare the models. In this case, the quadratic and cubic terms are not statistically significant themselves nor are their models statistically significant beyond the linear model. However, in a real research study, there would be other practical considerations to make before deciding on a final model.

More On Interactions, Polynomials, and HLR

Certainly, much more can be done with these topics than I have covered in my tutorials. What I have provided is a basic discussion with guided examples. The regression topics covered in these tutorials can be mixed and matched to create exceedingly complex models. For example, multiple interactions and higher order variables could be contained in a single model. The good news is that more complex models can be created using the same techniques covered here. The basic principles remain the same.

Complete Polynomial Regression Example

To see a complete example of how polynomial regression models can be created in R, please download the polynomial regression example (.txt) file.

R Tutorial Series: Regression With Categorical Variables

noreply@blogger.com (John Quick) — Mon, 01 Feb 2010 14:00:00 +0000

Categorical predictors can be incorporated into regression analysis, provided that they are properly prepared and interpreted. This tutorial will explore how categorical variables can be handled in R.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. Note that all code samples in this tutorial assume that this data has already been read into an R variable and has been attached. This dataset contains variables for the following information related to NFL quarterback and team salaries in 1991.

TEAM: Name of team
QB: Starting quarterback salary in thousands of dollars
TOTAL: team salary in thousands of dollars
CONF: conference (NFC or AFC)

In this dataset, the CONF variable is categorical. It can take on one of two values, either NFC or AFC. Suppose for the purposes of this tutorial that our research question is "how well do quarterback salary and conference predict total team salary?" The model that we use to answer this question will need to incorporate the categorical predictor for conference.

Dummy Coding

To be able to perform regression with a categorical variable, it must first be coded. Here, I will use the as.numeric(VAR) function, where VAR is the categorical variable, to dummy code the CONF predictor. As a result, CONF will represent NFC as 1 and AFC as 0. The sample code below demonstrates this process.

> #represent a categorical variable numerically using as.numeric(VAR)

> #dummy code the CONF variable into NFC = 1 and AFC = 0

> dCONF <- as.numeric(CONF) - 1

Note that the -1 that comes after the as.numeric(CONF) function causes the variables to read 1 and 0 rather than 2 and 1, which is the default behavior.

Interpretation

Visual

One useful way to visualize the relationship between a categorical and continuous variable is through a box plot. When dealing with categorical variables, R automatically creates such a graph via the plot() function (see Scatterplots). The CONF variable is graphically compared to TOTAL in the following sample code.

> #use the plot() function to create a box plot

> #what does the relationship between conference and team salary look like?

> plot(CONF, TOTAL, main="Team Salary by Conference", xlab="Conference", ylab="Salary ($1,000s)")

The resulting box plot is show below.

From a box plot, we can derive many useful insights, such as the minimum, maximum, and median values. Our box plot of total team salary on conference suggests that, compared to AFC teams, NFC teams have slightly higher salaries on average and the range of these salaries is larger.

Routine Analysis

Once a categorical variable has been quantified, it can be used in routine analyses, such as descriptive statistics and correlations. The following code depicts a few examples.

> #what are the mean and standard deviation of conference?

> mean(dCONF)

> [1] 0.5

> sd(dCONF)

> [1] 0.5091751

> #this makes sense… there are an even number of teams in both conferences and they are coded as either 0 or 1!

> #what is the correlation between total team salary and conference?

> cor(dCONF, TOTAL)

> [1]0.007019319

The correlation between total team salary and conference indicates that there is little to no linear relationship between the variables.

Linear Regression

Let's return to our original question of how well quarterback salary and conference predict team salary. With the categorical predictor quantified, we can create a regression model for this relationship, as demonstrated below.

> #create a linear model using lm(FORMULA, DATAVAR)

> #predict team salary using quarterback salary and conference

linearModel <- lm(TOTAL ~ QB + dCONF, datavar)

#generate model summary

summary(linearModel)

The model summary is pictured below.

Considering both the counterintuitive and statistically insignificant results of this model, our analysis of the conference variable would likely end or change directions at this point. However, there is one more interpretation method that is worth mentioning for future reference.

Split Model

With a dummy coded predictor, a regression model can be split into two halves by substituting in the possible values for the categorical variable. For example, we can think of our model as a regression of total salary on quarterback salary for two states of the world - teams in the AFC and teams in the NFC. These derivative models are covered in the following sample code.

> #input the categorical values to split the linear model into two representations

> #the original model: TOTAL = 19099 + 2.5 * QB - 103 * dCONF

> #substitute 0 for dCONF to derive the AFC model: TOTAL = 19099 + 2.5 * QB

> #substitute 1 for dCONF to derive the NFC model: TOTAL = 18996 + 2.5 * QB

#what is the predicted salary for a team with a quarterback salary of $2,000,000 in the AFC and NFC conferences?

#AFC prediction

19099 + 2.5 * 2000

[1] 24099

#NFC prediction

18996 + 2.5 * 2000

[1] 23996

Based only on what we have modeled, we can further infer that conference was not a significant predictor of total team salaries in the NFL in 1991. The difference between the team salaries based on conference is less than one-half of one percent on average! Of course, only using quarterback salary and conference to predict an NFL team's overall salary is neglecting quite a few potentially significant predictors. Nonetheless, split model interpretation is a useful way to break down the perspectives captured by a categorical regression model.

More On Categorical Predictors

Certainly, much more can be done with categorical variables than the basic dummy coding that was demonstrated here. Individuals whose work requires a deeper inspection into the procedures of categorical regression are encouraged to seek additional resources (and to consider writing a guest tutorial for this series).

Complete Categorical Regression Example

To see a complete example of how a categorical regression model can be created in R, please download the categorical regression example (.txt) file.

References

The Associated Press. (1991). Q-back and team salaries [Data File]. Retrieved December 14, 2009 from http://lib.stat.cmu.edu/DASL/Datafiles/qbacksalarydat.html