Receiver Operating Characteristic (ROC Curve) (SPSS)

What is a Receiver Operating Characteristic, or ROC Curve, as it is more commonly referred to as, and why should you care? The information provided within this article requires a firm understanding of the prior articles pertaining to logistic regression.

The ROC curve was first utilized during World War II for the analysis of radar signals.* Today, the ROC curve is utilized to illustrate, through the use of a diagram, the measurement of positive predictive results against false positive results. If this sounds complicated, don’t be nervous, the concept of the ROC curve is rather synthetic, meaning, that it is not at inherently intuitive. Therefore, I will do my best to explain its interpretation throughout this entry.

This article utilizes the prior the data set example and output.

To create an ROC curve, select “Analyze” from the drop down menu above the data sheet, then select “ROC Curve”.

Once this has been accomplished, you must specify your “Test Variable” and “State Variable” within the menu interface. In our case, the Test Variable will be: “Predicted Probability”, and the State Variable will be: “Cancer”. You must specify the positive value of the “State Variable”, meaning, that you must specify the value of the “State Variable” which identifies a positive result. This value will typically be “1”, which I have entered into the interface example below. For “Display” options, I have selected: “ROC Curve”, “With diagonal reference line”, “Coordinate points of the ROC Curve”.

This produces the output:

“Case Processing Summary” – This output demonstrates the number of actual positive and negative results contained within the original data set.

“ROC Curve” - What you are witnessing in the diagram above, is sensitivity vs. specificity as it pertains to points within the “Coordiantes of the Curve”. Y = Sensitivity and 

X = 1 – Specificity. I am aware that the graphical co-ordinates are plotted backwards on the graphic, however, I wanted to maintain the order in which they were presented within the table. What the diagram seeks to illustrate, is the predictive capacity of our model at various percentage confidence thresholds. Ideally, we would like to see a ROC curve illustration which contains single point in the upper leftmost corner of the graph. The green line illustrates a perfect binomial outcome (50/50 chance). If our blue line passes below the green line, the point in which this occurs, indicates a model cutoff which would provide probability results worse with than random chance. 

The ROC Curve graphic is providing an illustration of the "Coordinates of the Curve".

What is occurring in this table, is the measurement of positive results and false positive results at each percentage output provided by the logistic regression model.

So, to break this down in simpler terms, assuming that we utilized the value on the right: “Positive if Greater Than or Equal to”, as a cutoff in which to deem all cases of that value or higher as positive results (1), “Sensitivity” would then indicate the number of cases which were actually positive and identified as such (when 1 is predicted and the result is actually 1), “1-Specificity” would indicate the number of cases which were identified as positive by the model, but were actually negative (when 1 is predicted by the result is actually 0).

As an example: While identifying true outcomes based on a model specificity of .806, we will correctly identify 25% of positive outcomes. Therefore, 75% of the total positive outcomes will be overlooked. Additionally, 14.3% outcomes will be identified as positive while they are actually negative. That is to say, that this would be the result if we applied our model to our sample data and considered any result with a logit(p) <= .806 as being negative (0), and any result with a logit(p) >= as being positive (1).

With this in mind, it would not be an exaggeration to refer to “Sensitivity” as the number of positive cases identified by the model (as a percentage). Additionally, we could refer to “1-Specificity” as the number of false positive cases identified by the model (as a percentage).

“Area Under the Curve” – Commonly abbreviated as AUC, this value is representative of exactly what the name indicates. AUC, if we were shading in our graphic, would resemble the following:

For this reason, you can probably assume why you would prefer a higher AUC value, as opposed to the alternative. AUC values cannot, for obvious reasons, exceed the value of 1.

I hope that this entry de-mystified the concept of the ROC curve, and hopefully, has provided you with the confidence and information needed to implement its usage into your own work. Until next time, stay tuned, Data Heads!

* https://en.wikipedia.org/wiki/Receiver_operating_characteristic

Correlated Variables (SPSS)

In a prior article, I discussed correlation as it pertained to the “R” platform. In this entry, I will demonstrate how to generate a similar summary within the “SPSS” platform. We will be utilizing the data from the previous exercise for our example.

From the “Analyze” option on the top menu, select “Correlate” and then “Bivariate”.

Select the variable that will be analyzed. There is no need to select any additional options.

This produces the output below.

Typically, Pearson Correlation (“r”) values are interpreted as follows:

No Association: r value of 0

Small Association: r value of .1 to .3 (or) -.1 to .3

Medium Association: r value .3 to .5 (or) -.3 to -.5

Large Association: r value of .5 to 1 (or) -.5 to -1

Perfect Association: r value of 1

A rough rule of thumb is that only variables that possess a large association should be further examined within the context of the model.

That’s all for now, Data Heads! In the next article, we will discuss the ROC Curve and how it pertains to logistic regression.

(R) Logistic Regression Analysis (Non-Binary Categorical Variables) (SPSS)

In a previous article we covered how to analyze data through the utilization of the logistic regression model. In the example that was presented, categorical data was conveniently binary in every instance. In this example, I will demonstrate how to utilize the logistic regression model when categorical data contains multiple categories.

We will again refer to the example data set below. I have added an additional category that specifies the “Race” of each individual surveyed.

With category labels enabled, the data resembles:

Performing the analysis is similar to the prior example, except, in this case, we will be selecting the “Categorical” option. After doing such, we will specify that “Race” is a categorical covariate. 

Moving forward with our analysis, we will receive the following as a portion of our output:

“Race” has been split into 4 separate variables, with “Race” as a single variable, remaining for evaluation as whole. 

Race(1) refers to the “Race” category: “White”.

Race(2) refers to the “Race” category: “African American”.

Race(3) refers to the “Race” category: “Asian”.

Race(4) refers to the “Race” category: “Indian”.

The “Race” category “Native American” is still accounted for within the context of the model. However, its value is that of the constant in addition to all other variables.

In this example case, our equation would resemble:

Logit(p) = -2.13055 + (Age * 0.03335) + (Obese * -0.56859) + (Smoking * 3.02867) + (White * -1.10077) + (African_American * -1.05379) + (Asian * -1.22213) + (Indian * 0.69143)

So, if we wanted to test our model probability for an individual who was:

55 Years of Age

Obese

A Smoker

White

The equation would resemble:

Logit(p) = -2.13055 + (55 * 0.03335) + (1 * -0.56859) + (1 * 3.02867) + (1 * -1.10077) + (0 * -1.05379) + (0 * -1.22213) + (0 * 0.69143)

So our logit(p) value would be: 1.06301

Which equals a positive probability of: 0.7432653

Additionally, if our model was tested for an individual who was:

26 Years of Age

Not Obese

A Smoker

Native American

Our equation would resemble:

Logit(p) = -2.13055 + (26 * 0.03335) + (0 * -0.56859) + (1 * 3.02867) + (0 * -1.10077) + (0 * -1.05379) + (0 * -1.22213) + (0 * 0.69143) 

Logit(p) would equal: 1.76522

Which equals a positive probability of: 0.8538622

You can test this model in R with the following code:

# Model Test Code #

Age <- 0

Obese <- 0

Smoking <- 0

White <- 0

African_American <- 0

Asian <- 0

Indian <- 0

p <- -2.13055 + (Age * 0.03335) + (Obese * -0.56859) + (Smoking * 3.02867) + (White * -1.10077) + (African_American * -1.05379) + (Asian * -1.22213) + (Indian * 0.69143)

plogis(p) 

Here is how you would create the same model through the utilization of the “R” Platform:

# Non-Binary Categorical Variables #

Age <- c(55, 45, 33, 22, 34, 56, 78, 47, 38, 68, 49, 34, 28, 61, 26)

Obese <- c(1,0,0,0,1,1,0,1,1,0,1,1,0,1,0)

Smoking <- c(1,0,0,1,1,1,0,0,1,0,0,1,0,1,1)

Cancer <- c(1,0,0,1,0,1,0,0,1,1,0,1,1,1,0)

White <- c(1,1,1,0,0,0,0,0,0,0,0,0,0,0,0)

African_American <- c(0,0,0,1,1,1,0,0,0,0,0,0,0,0,0)

Asian <- c(0,0,0,0,0,0,1,1,1,0,0,0,0,0,0)

Indian <- c(0,0,0,0,0,0,0,0,0,1,1,1,0,0,0)

Native_American <- c(0,0,0,0,0,0,0,0,0,0,0,0,1,1,1)

CancerModelII <- data.frame(Age, Obese, Smoking, Cancer, White, African_American, Asian, Indian, Native_American )

CancerModelLogII <- glm(Cancer~ Age + Obese + Smoking + White + African_American + Asian + Indian + Native_American, family=binomial)

summary(CancerModelLogII)

# Which produces the output #

Call:
glm(formula = Cancer ~ Age + Obese + Smoking + White + African_American +
Asian + Indian + Native_American, family = binomial)

Deviance Residuals:
      Min       1Q      Median      3Q        Max
-1.9613 -0.7252    0.4240    0.8107   1.7092

Coefficients: (1 not defined because of singularities)
                                             Estimate        Std. Error      z value     Pr(>|z|)
(Intercept)                           -2.13055         2.58207         -0.825       0.409
Age                                      0.03335         0.04641           0.719       0.472
Obese                                 -0.56859         1.60680          -0.354       0.723
Smoking                              3.02867          1.95858          1.546        0.122
White                                 -1.10077          2.35673          -0.467       0.640
African_American             -1.05379          2.18843          -0.482       0.630
Asian                                 -1.22213          2.40838          -0.507       0.612
Indian                                 0.69143          2.51153            0.275       0.783
Native_American                 NA                 NA                  NA             NA

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 20.728 on 14 degrees of freedom
Residual deviance: 15.366 on 7 degrees of freedom
AIC: 31.366

Number of Fisher Scoring iterations: 4

(R) Logistic Regression Analysis (Binary Categorical Variables) (SPSS)

Today we will be discussing an advanced topic, but a useful topic nonetheless, that topic being: Logistic Regression. Before delving into this subject matter, I would advise you, the reader, if you are not already familiar with linear regression, to please review the articles pertaining to this topic.

Logistic regression is a method which is similar to linear regression, however, the logistic regression method is utilized specifically to create models which analyze dependent variables that are binary in nature (True/False, Yes/No, Positive/Negative, 1/0, etc.).

Let’s begin with a simple example, I have previously entered some sample data into SPSS.

As we discussed in a previous article, you should edit the variable labels so that they correspond with the appropriate binary outcomes.

To perform a logistic regression analysis, select "Analyze" from the top drop down menu, then select "Regression". From the next menu, select "Binary Logistic".

The next screen presents us with the options necessary to structure the equation utilized to create the model. The logistic regression model is structured in a way that is similar to the linear regression model. In the case of our example, “Cancer” will be our dependent variable. “Age”, “Obese”, and “Smoking” will be the model’s independent variables.

NOTE: In this example, our categorical variables are binary. Meaning, that they can only represents two values (0,1). If we were working with sample data that contained categorical variables that were not binary, we would have to specify this in the following screen. This menu can be toggled by clicking on the “Categorical” option.

From the options menu, select “Hosmer-Lemeshow goodness-of-fit”. This test will be useful in interpretating the model results.

Selecting “Save” from the menu options, will present you with this interface. Check “Probabilities” and “Group Membership” beneath the “Predicted Values“ menu header. These options being enabled will output probability data directly to our SPSS spreadsheet.

So what does all of the generated output indicate?

For the most part, you can ignore all of the data listed in the screen below:

Block 0

(Ignore Block 0)

Block 1

Starting from Block 1 to perform analysis, please first observe the portion of this output which reads “Omnibus Tests of Model Coefficients”. We will primarily be concerned with the bottom-most line which coincides with the row labeled, "Model".

The rightmost column entry: "Sig.", indicates the significance level of the model. Depending on the pre-determined confidence interval which has been established, this value will determine whether or not the model is significant. "df" represents the test model's degrees of freedom. The "Chi-square" column contains the corresponding test statistic. All three values, in addition to the number of observations within the model, are utilized in tandem when documenting the model's significance.    

“-2 Log Likelihood”, “Cox & Snell R Square”, “Nagelkerke R Square” – All three columns correspond with outputs generated from methodologies which test the strength of the model.

Cox & Snell R Squared – This value can reach a maximum of .75. Therefore, I believe that it is best to consult the “Nagelkerke R Square” value when assessing the overall model.

Nagelkerke R Square – This value can reach a maximum of 1. Thus, it is the measurement that best resembles the equivalent of The Coefficient of Determination. I would recommend referring to this value when considering the strength of the predictive model.

The “Homer and Lameshow Test” – This is measuring for co-linearity amongst the independent variables of the model. Only the significance of this test is valuable. If the significance percentage is lower than .05, a correlation may be present. In such cases,  additional testing should be performed to address such. Since this test utilizes a Chi-Squared distribution, it is sensitive to sample size. The test is considered most accurate when a sample size is greater than 400.

“Contingency Table for Hosmer and Lameshow Test” – This can be ignored.

“Classification Table” – This output is simple to understand but difficult to follow. What is being illustrated by the output is the number of positive and negative cases predicted by the model, and the number of actual cases that occurred within the sample.

For our model, there were 5 predicted non-occurrences of cancer, and 2 occurrence of cancer that were not predicted. Thus, the percentage of correct estimates is 71.4. (5/(5+2))

Also, there were 2 predicted occurrences of cancer, and 6 non predicted occurrence. Thus, the percentage of correct estimates is 75. (6/(6+2))

The total percentage of accurate estimations provided by the model is 73.3. ((6+5)/(6+2+5+2).

“Variables in the Equation”

The first column lists the model variables. “Age”, “Obese”, “Smoking” are dependent variables. “Constant” is a variable which will be included as an aspect of the model.

So, if we were to construct the model as an equation, it would resemble:

Logit(p) = (Age * .030) + (Obese * -.389) + (Smoking * 2.544) – 2.344

S.E - represents the standard error of each variable within the context of the model.

Wald - represents the value of the Wald Test variable. You can ignore this value.

df - is indicating the degrees of freedom for each variable.

Sig -  is indicating the significance of each variable within the context of the model. An insignificant variable value (typically > .05), does not necessarily indicate that the variable itself should be excluded from subsequent model generation.

Exp(B) or “Odds Ratio” – This indicates the value in which the probability of a positive outcome will increase if the corresponding variable is increased by one. This warrants further explanation, as logistic regression models utilize the “Log” function to generate results.

Returning to our original model equation:

Logit(p) = (Age * .030) + (Obese * -.389) + (Smoking * 2.544) – 2.344

Let’s assume that we wish to test our model on a 19 year old, who is not obese, and who smokes. The equation would resemble:

Logit(p) = (19 * .030) + (0 * -.389) + (1 * 2.544) – 2.344

Logit(p) = .77

But Logit(p) does not equal probability. It equals the Logit value of .77.

To generate the associated probability of this value, you can enter the following information into the R console:

plogis(.77)

Which produces the console output:

[1] 0.6835209

Additionally, you could also create the model within “R” by utilizing the code:

a <- 19
b <- 0
c <- 1

p <- (.030 * a) + (-.389 * b) + (2.544 * c) - 2.344

plogis(p)

Which produces the console output:

[1] 0.6835209

In both cases, the output value is the probability value of a positive occurrence. In the case of our scenario, this would equate a cancer diagnosis. The numbers from our model are un-realistic, but this only due to the absolutely un-realistic sample data that I hastily created.

Returning to the odds ratio.

The value given in the column Exp(B) is the exponential value of the value given in column B. What this value is indicates is the probability increase, within the context of the model, if the variable value of the corresponding B variable is increased by one. However, this value must first be transformed.

Assuming the context of the model, if an individual were being assessed for cancer, assuming that he aged a year, the probability of cancer would increase:

(1.031 – 1) * 100

3.1 %

Data Sheet Output

Previously, during this exercise, I asked that you select “Save” from the menu options, and then subsequently select “Probabilities” and “Group Membership”. In doing so, we have informed SPSS to output data directly to our datasheet. The results are listed below.

We are presented with two new variables “PRE_1” and “PGR_1”. “PRE_1” represents the probability (.00 – 1.00) of a positive event occurring when the dependent variable data found in the left most columns is input into the model. “PGR_1” represents the model’s prediction given the dependent variable data (1 or 0).  This data output is needed to create a ROC curve. ROC curves will be covered in a later article.

Model Creation within R

I will now briefly illustrate how to obtain the same results in R.

# Model Creation #

Age <- c(55, 45, 33, 22, 34, 56, 78, 47, 38, 68, 49, 34, 28, 61, 26)

Obese <- c(1,0,0,0,1,1,0,1,1,0,1,1,0,1,0)

Smoking <- c(1,0,0,1,1,1,0,0,1,0,0,1,0,1,1)

Cancer <- c(1,0,0,1,0,1,0,0,1,1,0,1,1,1,0)

CancerModel <- data.frame(Age, Obese, Smoking, Cancer) 

# Analyze the Significance of the Model #

CancerModelLog <- glm(Cancer ~ Age + Obese + Smoking, family=binomial)

CancerModelLog1 <- glm(Cancer ~ 1, family=binomial)

anova(CancerModelLog, CancerModelLog1, test="LRT")

# Console Output #

Analysis of Deviance Table

Model 1: Cancer ~ Age + Obese + Smoking
Model 2: Cancer ~ 1
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 11 16.807
2 14 20.728 -3 -3.9209 0.2701

# Summary Creation and Output # 

CancerModelLog <- glm(Cancer~ Age + Obese + Smoking, family=binomial)

summary(CancerModelLog)

# Output #

Call: 

glm(formula = Cancer ~ Age + Obese + Smoking, family = binomial)

Deviance Residuals:
    Min         1Q          Median     3Q        Max
-1.6096     -0.7471     0.5980    0.8260  1.8485

Coefficients:
                      Estimate      Std. Error    z value     Pr(>|z|)
(Intercept)     -2.34431       2.25748        -1.038      0.2991
Age                0.02984       0.04055          0.736      0.4617
Obese            -0.38924      1.39132         -0.280      0.7797
Smoking        2.54387       1.53564          1.657      0.0976 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 20.728 on 14 degrees of freedom
Residual deviance: 16.807 on 11 degrees of freedom
AIC: 24.807
Number of Fisher Scoring iterations: 4

# Generate Nagelkerke R Squared #

# Download and Enable Package: "BaylorEdPsych" #

PseudoR2(CancerModelLog)

# Console Output # McFadden   Adj.McFadden    Cox.Snell    Nagelkerke    McKelvey.Zavoina      Effron 0.2328838    -0.2495624         0.2751639    0.3674311    0.3477522               0.3042371 0.8000000 Adj.Count          AIC          Corrected.AIC 0.5714286      23.9005542      27.9005542 # Calculate Exp(B) # # (Intercept) # exp(-2.34431) # Age # exp(0.02984) # Obese # exp(-0.38924) # Smoking # exp(2.54387) # Output # [1] 0.09591336 [1] 1.03029 [1] 0.6775716 [1] 12.72884 # Utilize VIF() and COR() in lieu of The Homer and Lameshow Test #

# Generate Correlation Matrix # 

correlationmatrix <- cor(CancerModel)

# Output #

                              Age             Obese          Smoking       Cancer
Age                 1.00000000     0.1231756  -0.2836428    0.02147097
Obese             0.12317556    1.0000000   0.4642857    0.19642857
Smoking       -0.28364280     0.4642857  1.0000000     0.46428571
Cancer          0.02147097    0.1964286   0.4642857    1.00000000

# Generate Variance Influence Factor #

# Download and Enable Package: "car" #

vif(CancerModelLog)

# Output # 

     Age        Obese    Smoking
1.286936 1.361140 1.654232

(R) Post Hoc Analysis and One Way ANOVA (SPSS)

As was previously mentioned, new entries posted on this blog will be primarily non-R related.

Today’s post will discuss Post Hoc Analysis, specifically Tukey’s Honest Significance Test. This test is also known as The Tukey Method, Tukey’s HSD, or TukeyHSD() in R.

Post Hoc refers to the testing that is performed following an ANOVA test. What this testing seeks to discover, is the significance of relationships that exist between variables within an ANOVA model. There are many different Post Hoc tests that can be utilized. For the purpose of this article, we will be specifically discussing Tukey’s HSD.

Something that I should mention before proceeding, is the reason for the utilization of ANOVA as opposed to a T-Test. ANOVA allows us to compare the means between various groups simultaneously, while maintaining the same confidence interval. If we had four experimental groups to test between, this would require 6 T-Tests.

1 vs. 2 | 1 vs. 3 | 1 vs. 4

2 vs. 3 | 2 vs. 4

3 vs. 4

Each T-Test, if assuming an alpha of .05, has a 5% chance of a Type I error occurring. This means, that there is a 30% chance (.05 * 6), that at least one Type I error would occur. The T-Test will analyze for a statistical difference between the means of two groups, whereas the ANOVA, analyzes for differences within the set of means.

If you recall from the previous article, we addressed two separate scenarios, one in which a cook was testing for the salt content of soup, and the other, in which the impact of study time was being assessed as it applied to students from two different schools.

We will run a Tukey’s HSD on the data collected from each study.

Scenario A: The Soup Scenario

satisfaction <- c(4, 1, 8, 4, 5, 3, 5, 3, 2, 5)

salt <- c(rep("low",3), rep("med",4), rep("high",3))

salttest <- data.frame(satisfaction, salt)

results <- aov(satisfaction~salt, data=salttest)

Now to run the Tukey HSD Post Hoc Inquiry:

TukeyHSD(results)

Which produces the output:

Tukey multiple comparisons of means

    95% family-wise confidence level

Fit: aov(formula = satisfaction ~ salt, data = salttest)

$salt

                          diff                lwr      upr          p adj

low-high  1.00000000 -4.148005 6.148005 0.8387911

med-high  0.91666667 -3.898852 5.732185 0.8445186

med-low  -0.08333333 -4.898852 4.732185 0.9985693

Let us review each aspect of this output:

Diff – Is the difference of the averages between the values.

Lwr – Is the lower confidence interval of the difference.

Upr – Is the upper confidence interval of the difference.

P adj – The p-values pertaining to the significance of the compound values. Again, if 95%, we will be looking for values of significance that are less than .05.

Each value within the “p adj” column corresponds to an assessment of the significance pertaining to separate categorical aspects of the model. In the case of the above output, assuming an alpha value of .05, there were no significant differences between any of the categorical factors (p = 0.839; p = 0.844; p = 0.999).

Scenario B: Schools, Study Time and Stress Scenario

satisfaction <- c(7, 2, 10, 2, 2, 8, 5, 1, 3, 10, 9, 10, 3, 10, 8, 7, 5, 6, 4, 10, 3, 6, 4, 7, 1, 5, 5, 2, 2, 2)

studytime <- c(rep("One Hour",10), rep("Two Hours",10), rep("Three Hours",10))

school = c(rep("SchoolA",5), rep("SchoolB",5), rep("SchoolA",5), rep("SchoolB",5), rep("SchoolA",5), rep("SchoolB",5))

schooltest <- data.frame(satisfaction, studytime, school)

results <- aov(lm(satisfaction ~ studytime * school, data=schooltest))

summary(results)

Now to run the Tukey HSD Post Hoc Inquiry:

TukeyHSD(results)

This produces the following output:

$studytime

                                          diff        lwr             upr      p adj

Three Hours-One Hour  -1.3 -4.5013364 1.901336 0.5753377

Two Hours-One Hour     2.2 -1.0013364 5.401336 0.2198626

Two Hours-Three Hours  3.5  0.2986636 6.701336 0.0302463

Is describing the relationship between the varying levels of study time as it pertains to stress.

The next portion of the output:

$school

                                diff       lwr             upr      p adj

SchoolB-SchoolA -0.6 -2.760257 1.560257 0.571817

Describes the relationship between the two school types as it pertains to stress.

Finally, the last portion of the output:

$`studytime:school`

                                                                        diff         lwr       upr             p adj

Three Hours:SchoolA-One Hour:SchoolA    -0.4  -6.005413 5.2054132 0.9999178

Two Hours:SchoolA-One Hour:SchoolA       3.4  -2.205413 9.0054132 0.4401459

One Hour:SchoolB-One Hour:SchoolA        0.8  -4.805413 6.4054132 0.9976117

Three Hours:SchoolB-One Hour:SchoolA    -1.4  -7.005413 4.2054132 0.9696463

Two Hours:SchoolB-One Hour:SchoolA       1.8  -3.805413 7.4054132 0.9157375

Two Hours:SchoolA-Three Hours:SchoolA    3.8  -1.805413 9.4054132 0.3223867

One Hour:SchoolB-Three Hours:SchoolA     1.2  -4.405413 6.8054132 0.9844928

Three Hours:SchoolB-Three Hours:SchoolA -1.0  -6.605413 4.6054132 0.9932117

Two Hours:SchoolB-Three Hours:SchoolA    2.2  -3.405413 7.8054132 0.8260605

One Hour:SchoolB-Two Hours:SchoolA      -2.6  -8.205413 3.0054132 0.7067715

Three Hours:SchoolB-Two Hours:SchoolA   -4.8 -10.405413 0.8054132 0.1240592

Two Hours:SchoolB-Two Hours:SchoolA     -1.6  -7.205413 4.0054132 0.9470847

Three Hours:SchoolB-One Hour:SchoolB    -2.2  -7.805413 3.4054132 0.8260605

Two Hours:SchoolB-One Hour:SchoolB       1.0  -4.605413 6.6054132 0.9932117

Two Hours:SchoolB-Three Hours:SchoolB    3.2  -2.405413 8.8054132 0.5052080

Describes the relationships between the combination of hours studied and school types.

We can make the following interpretations from the above outputs:

There was a significant difference in stress level between students who study two hours and students who study three hours (p = 0.0302463).

There was not a significant difference in stress level between students who attend SchoolA, and students who attend SchoolB.

There were not a significant differences in stress levels as it pertains to the combination of factors: school and study time.

I have often been asked what differentiates an ANOVA Post Hoc Test (such as Tukey’s HSD), from a T-Test proceeding an ANOVA calculation. The reasons for performing a Post Hoc Test, Tukey's in our case, as opposed to a T-Test, are as follows:

1. Performing multiple T-Tests to check for significance is ultimately time consuming, and nullifies the initial convenience of running an ANOVA test. Additionally, doing such, re-creates the compounding probability of error that we originally sought to avoid.

2. Tukey’s HSD takes into account the significance of each variable as they interact with other variables within the ANOVA model. Re-testing with the T-Test may show what data sets INDENPENDENTLY differ from the other data sets, but it will not illustrate what data sets differ within the model.

Now, I will demonstrate how to perform both a One Way ANOVA Test and a Post Hoc Tukey’s HSD test within SPSS .

First, we will need to define the variables, this can be achieved within the “Variable View” portion of SPSS. Selecting this view can be achieved by clicking the “Variable View” tab on the lower right hand side of the SPSS console.

Once “Variable View” has been selected, we can begin by defining our variable types.



Here I have defined two variables, “Satisfaction” and “Soup”, both were assigned the default variable parameters by the SPSS system.

Next, we need to define our value labels, to achieve this, I clicked on the cell which coincides with the variable “Soup”. This brings up a user interface, which allows for the entry of value labels, and the value for which the label is assigned.



Once this data has been input, we can now input the corresponding values into SPSS which are required for the assembly of our ANOVA model.

When this step has been completed, to proceed, we must choose “Analyze” from the upper most drop down menu. Select the option “Compare Means”, and the subsequent option, “One-Way ANOVA”.

This course of action should cause a menu to appear. For our “Dependent List” variable, we will choose, “Satisfaction”. For our “Factor” variable, we will choose “Soup”. Once this has been completed, select the middle box on the right corner of the menu which reads, “Post Hoc”. This causes another menu to appear which presents Post Hoc Analysis options. For our purposes, we will be checking the box next to “Tukey” prior to proceeding. Significance level should be left at .05 (or Alpha = .05). Click “Continue”, then click, “OK”.

This presents a more detailed Tukey’s HSD output than what was originally available in R:

Compared to the R output, the following output is detailing:

Mean Difference – Is the difference of the averages between the values. “diff” in R.

Std. Error – The standard error of the compared values. No equivalent in R.

Lower Bound – Is the lower confidence interval of the difference. “lwr” in R.

Upper Bound – Is the upper confidence interval of the difference. “upr” in R.

Sig. – The p-values pertaining to the significance of the compared values. Again, if 95%, we will be looking for values of significance that are less than .05. “p adj” in R.

That’s all for now, Data Heads. I’ll see you again soon with a brand new article, the subject matter of such is undetermined.

(R) Analysis of Variance - ANOVA

In this article, we will discuss ANOVA, specifically, when its usage is appropriate, and how it can be utilized within R. This will likely be the final article of the current series of entries pertaining to The R Programming Language. Subsequent articles will discuss concepts and usage of software within the SPSS platform.

ANOVA is an abbreviation that represents a method known as The Analysis of Variance.

There are few terms that are specific to ANOVA, those are:

Way – Which refers to an independent variable within the ANOVA model.

Factor – Another term which refers to an independent variable.

Level – The category of an independent variable within the ANOVA model.

ANOVA is used to compare the variances of various sample groups against one another. In many ways it is similar to a t-test, however, ANOVA allows for multiple group comparisons. This differs from the t-test, which only allows for one single group to be compared to another single group.

A post-hoc test is often performed after ANOVA has been calculated. We will discuss this topic in a different article. A post-hoc test is used to further investigate data sample similarities and is utilized when the ANOVA model returns certain results.

Like the t-test, there are different variations of the ANOVA model that are applicable depending on the data being analyzed. We will review three common ANOVA application as they pertain to various data types. The analyzation of the output of the model data is performed through the utilization of the F-Test. For a detailed description of the F-Test, and what conclusions it provides, please refer to the pervious article.

One Way ANOVA 

As a reminder, Way, in this scenario, is referring to a single independent variable.

In a one way ANOVA, we are assuming the following:

1. Each sample is random.
2. Each sample is in no way influenced by the other sampling results.
3. Each dependent variable is sampled from a normally distributed population.
4. The variances of the samples, should be equivalent, or somewhat equivalent. The reason for such, is that the population variances are assumed to be equal for each sample.

The hypothesis for this model type will be:

H0: u1 = u2 = u3 =…..etc.

H1: Not all means are equal.

Example Problem:

A chef wants to test if patrons prefer a soup which he prepares based on salt content. He prepares a limited experiment in which he creates three types of soup: soup with a low amount of salt, soup with a high amount of salt, and soup with a medium amount of salt. He then servers this soup to his customers and asks them to rate their satisfaction on a scale from 1-8.

Low Salt Soup it rated: 4, 1, 8
Medium Salt Soup is rated: 4, 5, 3, 5
High Salt Soup is rated: 3, 2, 5

Hypothesis:

H0: u1 = u2 = u3 =…..etc.

H1: Not all means are equal.

Let’s use this data to create a model within R:

satisfaction <- c(4, 1, 8, 4, 5, 3, 5, 3, 2, 5)

salt <- c(rep("low",3), rep("med",4), rep("high",3))

salttest <- data.frame(satisfaction, salt)

results <- aov(satisfaction~salt, data=salttest)

summary(results)

This produces the output:

                 Df     Sum Sq Mean Sq     F value   Pr(>F)
salt             2       1.92     0.958           0.209     0.816
Residuals   7       32.08   4.583   

If p < .05, we will reject the null hypothesis.

Hypothesis: 0.816 > .05

Since the model’s p-value (.816) is greater than the assumed alpha (.05), we will fail to reject the null hypothesis. What this is indicating, is that at 95% confidence interval, we cannot state that through the analysis of the data provided, that there is a significant difference of customer satisfaction as it pertains to salt content in soup.

Two Way ANOVA

Two way, in this scenario, is referring to the two independent variables which will be utilized within this ANOVA model.

The hypothesis for this model type will be:

1.

H0: u1 = u2 = u3 =…..etc. (All means are equal)

H1: Not all means are equal.

2.

H0: uVar1 = uVar2 (Var1’s value does not significantly differ from Var2’s value)

H1: uVar1 NE uVar2

3.

H0: An interaction is absent.

H1: An interaction is present.

Example Problem:

Researchers want to test study habits within two schools as they pertain to student life satisfaction. The researchers also believe that the school that each group of students is attending may also have an impact on study habits. Students from each school are assigned study material which in sum, totals to 1 hour, 2 hours, and 3 hours on a daily basis. Measured is the satisfaction of each student group on a scale from 1-10 after a 1 month duration.

School A:

1 Hour of Study Time: 7, 2, 10, 2, 2
2 Hours of Study Time: 9, 10, 3, 10, 8
3 Hours of Study Time: 3, 6, 4, 7, 1

School B:

1 Hour of Study Time: 8, 5, 1, 3, 10
2 Hours of Study Time: 7, 5, 6, 4, 10
3 Hours of Study Time: 5, 5, 2, 2, 2

Let’s state our hypothesizes, as they apply to this problem:

1.

H0: u1 = u2 = u3 (Stress levels DO NOT differ depending on hours of daily study.)

H1: Not all means are equal. (Stress levels DO differ depending on hours of daily study.)

2.

H0: uSchoolA = uSchoolB (Stress levels DO NOT significantly differ depending on school school.)

H1: uSchoolA NE uSchoolB (Stress levels DO significantly differ depending of school.)

3.

H0: An interaction is absent. (The combination of school and study time is NOT impacting the outcome)

H1: An interaction is present. (The combination of school and study time IS impacting the outcome)

Entering this into R can be tricky, but stay with me:

satisfaction <- c(7, 2, 10, 2, 2, 8, 5, 1, 3, 10, 9, 10, 3, 10, 8, 7, 5, 6, 4, 10, 3, 6, 4, 7, 1, 5, 5, 2, 2, 2)

studytime <- c(rep("One Hour",10), rep("Two Hours",10), rep("Three Hours",10))

school = c(rep("SchoolA",5), rep("SchoolB",5), rep("SchoolA",5), rep("SchoolB",5), rep("SchoolA",5), rep("SchoolB",5))

schooltest <- data.frame(satisfaction, studytime, school)

results <- aov(lm(satisfaction ~ studytime * school, data=schooltest))

summary(results)

Which produces the output:

                                 Df      Sum Sq      Mean Sq      F value     Pr(>F)
studytime                  2          62.6          31.300        3.809        0.0366 *
school                      1            2.7            2.700         0.329        0.5718
studytime:school     2            7.8            3.900         0.475        0.6278
Residuals               24        197.2           8.217
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Since we have three hypothesis tests, we must assess all three of the p-values present within the output.

Study Time

p = 0.0366

School

p = 0.5718

Study Time : School

p = 0.6278

In investigating the output we can make the following conclusions:

Hypothesis 1: 0.0366 < .05

Hypothesis 2: 0.5718 > .05

Hypothesis 3: 0.6278 > .05

If p < .05, we will reject the null hypothesis.

Hypothesis 1: Reject

Hypothesis 2: Fail to Reject

Hypothesis 3: Fail to Reject

So we can state:

Students of different schools did not significantly different stress levels. There was significant difference between the levels of study time as it pertains to stress. No interaction effect was present.

(Two Way ANOVA must have columns observations of equal length)

Repeated-Measures ANOVA

A repeated measures ANOVA is similar to a paired t-test in that it samples from the same set more than once. This model contains one factor with at least two levels, and the levels are dependent.

Example Problem:

Researchers want to test the impact of reading existential philosophy on a group of 8 individuals. They measure the happiness of the participants three times, once prior to reading, once after reading the materials for one week, and once after reading the materials for two weeks. We will assume an alpha of .05.

Before Reading = 1, 8, 2, 4, 4, 10, 2, 9
After Reading = 4, 2, 5, 4, 3, 4, 2, 1
After Reading (wk. 2) = 5, 10, 1, 1, 4, 6, 1, 8

Hypothesis:

H0: u1 = u2 = u3

H1: Not all means are equal.

Let’s use this data to create a model within R:

library(lme4) # You will need to install and enable this package #

happiness <- c(1, 8, 2, 4, 4, 10, 2, 9, 4, 2, 5, 4, 3, 4, 2, 1, 5, 10, 1, 1, 4, 6, 1, 8 )

week <- c(rep("Before", 8), rep("Week1", 8), rep("Week2", 8))

id <- c(1,2,3,4,5,6,7, 8)

survey <- data.frame(id, happiness, week)

model <- lmer(happiness ~ week + (1|id), data=survey)

anova(model)

Which produces the output:

         Analysis of Variance Table
            Df      Sum Sq      Mean Sq      F value
week     2      15.083         7.5417        1.0462

The F-Test statistic = 1.0462
 
To calculate the p-value of our test statistic, we can use the following r-code:

pf(q=1.0462, df1=2, df2=14, lower.tail=FALSE) # Test Statistic , Numerator Degrees of Freedom = 2, Denominator Degrees of Freedom = 14 #

Which produces the output:

[1] 0.3771816

If p < .05, we will reject the null hypothesis.

Hypothesis: 0.3771816 > .05

With this information, we can conclude that the three conditions did not significantly differ pertaining to level of happiness. 

A similar methodology that can be utilized to perform this analysis:

library(lme4) # You will need to install and enable this package #
library(nlme) # You will also need to install and enable this package #

happiness <- c(1, 8, 2, 4, 4, 10, 2, 9, 4, 2, 5, 4, 3, 4, 2, 1, 5, 10, 1, 1, 4, 6, 1, 8 )

week <- c(rep("Before", 8), rep("Week1", 8), rep("Week2", 8))

id <- c(1,2,3,4,5,6,7, 8)

survey <- data.frame(id, happiness, week)

model <- lme(happiness ~ week, random=~1|id, data=survey)

anova(model)

This method saves some time by producing the output:

                   numDF   denDF    F-value         p-value
(Intercept)     1            14         37.21053        <.0001
week             2            14           1.04624          0.3772

That is all for now, Data Heads. The topic of the next article will be Post-Hoc Analysis. Stay tuned!

(R) F-Test

You may remember the F-Test from the previous article on multiple linear regression. In this entry, we will further delve into the concept of the F-Test.

The F-Test is a statistical method for comparing two population variances. It’s most recognized utilization is as one of the aspects of the ANOVA method. This method will be discussed in a later article.

Essentially, the F-Test model enables the creation of a test statistic, critical value and a distribution model. With these values derived, a hypothesis test can be stated, and from such, the comparison of two variance can be achieved.

Some things to keep in mind before moving forward:

1. The F-Test assumes that the samples provided, originated from a normal distribution.

2. The F-Test attempts to discover whether two samples originate from populations with equal variances.

So for example, if we were comparing the following two samples:

samp1 <- c(-0.73544189, 0.36905647, 0.69982679, -0.91131589, -1.84019291, -1.02226811, -1.85088278, 2.24406451, 0.63377787, -0.80777949, 0.60145711, 0.43853971, -1.76386879, 0.32665597, 0.32333871, 0.90197004, 0.29803556, 0.47333427, 0.23710263, -1.48582332, -0.45548478, 0.36490345, -0.08390139, -0.46540965, -1.66657385)

samp2 <- c(0.67033912, -1.23197505, -0.18679478, 1.06563032, 0.08998155, 0.22634414, 0.06541938, -0.22454059, -1.00731073, -1.43042950, -0.62312404, -0.22700636, -0.71908729, -0.36873910, 0.15653935, -0.19328338, 0.56259671, 0.31443699, 1.02898245, 1.18903593, -0.14576090, 0.68375259, -0.15348007, 1.58654607, 0.01616986)

For a right tailed test, we would state the following hypothesis:

H0: σ2/1 =σ2/2
Ha: σ2/1>σ2/2

# Null Hypothesis = Variances are equal. #

# Alternative Hypothesis = The first measurement of variance is greater than the second measurement of variance. #

With both samples imported into R, we can now utilize the following code to perform the F-Test:

(We will assume an alpha of .05):

var.test(samp1, samp2, alternative = "greater", conf.level = .95)

Which produces the output:

    F test to compare two variances

data: samp1 and samp2
F = 1.9112, num df = 24, denom df = 24, p-value =
0.05975
alternative hypothesis: true ratio of variances is greater than 1
95 percent confidence interval:
0.9634237 Inf
sample estimates:
ratio of variances 

        1.911201

Let us review each aspect of this output:

“ F = “ is the F-Test test statistic.

“num df = “ is the value of the degrees of freedom found within the numerator.

“denom df = “ is the value of the degrees of freedom found within the denomenator.

“p-value = “ is the probability of the corresponding F-Test statistic.

“95 percent confidence interval:” is the ratio between the two population variances at the 95% confidence level.

“ratio of variances” is the value of the variance of sample 1 divided by the variance of sample 2.

Looking at the p-value, which is greater than our alpha value (0.05975 > .05), we cannot conclude, that at a 95% confidence level, that our samples were taken from populations with differing variances.

Additionally, we can confirm this conclusions by comparing our F-Test statistic of 1.9112, to the F-Value which coincides with the appropriate degrees of freedom and alpha value. To find this value, we would typically consult a chart in the back of a statistics textbook. However, R makes the situation simpler by providing us with a method to reference this value.

Utilizing the code:

qf(.95, df1=24, df2=24) #Alpha .05, Numerator Degrees of Freedom = 24, Denomenator Degrees of Freedom = 24#

Produces the output:

[1] 1.98376

Again, we cannot conclude that because 1.9112 < 1.98376, that our samples were taken from populations with differing variances.

If we were to graph this test and distribution, the illustration would resemble:

If you would like to create your own f-distribution graphs, sans the mark-ups, you could use the following code:

curve(df(x, df1=24, df2=24), from=0, to=5) # Modify the degrees of freedom only #

Below is an illustration of a few various f-distribution types by varying degrees of freedom:

I hope that you found this article useful, in the next post, we will begin to discuss the concept of ANOVA.

* A helpful article pertaining to the F-Test statistic: http://atomic.phys.uni-sofia.bg/local/nist-e-handbook/e-handbook/eda/section3/eda359.htm

** Source for F-Distribution Image: https://en.wikipedia.org/wiki/F-distribution