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:
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
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:
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.
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.
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).
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.
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:
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.
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 #
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:
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.
Revisiting multiple regression for the purpose of establishing a greater understanding of the topic, and additionally, to foster a more efficient application of the regression model; in this entry we will review functions and statistical concepts which were overlooked in the prior article pertaining to the subject matter.
Variance Influence Factor or VIF():
What is VIF()? According to Wikipedia, it is defined as such: “The variance inflation factor (VIF) quantifies the severity of multicollinearity in an ordinary least squares regression analysis.”
To define VIF() in layman’s terms, The Variance Influence Factor is a method for weighing each variable’s coefficient of determination (R-Squared), against all other variables within a multiple regression equation. It’s really as simple as that, but none of this will make sense without an example.
Example:
Consider the following numerical sets: w <- c(23, 42, 55, 16, 24, 27, 24, 15, 23, 85) x <- c(27, 34, 22, 30, 17, 32, 25, 34, 46, 37) y <- c(70, 80, 73, 77, 60, 93, 85, 72, 90, 85) z <- c(13, 22, 18, 30, 15, 17, 20, 11, 20, 25)
Let’s utilize these values to create a multiple regression model:
lm.multiregressw <- (lm(w ~ x + y + z))
Now let’s take a look at that model:
summary(lm.multiregressw) Call: lm(formula = w ~ x + y + z)
Residuals: Min 1Q Median 3Q Max -29.701 -11.020 -4.462 3.465 44.108
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.40034 68.57307 -0.020 0.984 x -0.08981 1.41549 -0.063 0.951 y 0.19770 1.20528 0.164 0.875 z 1.15241 1.60560 0.718 0.500
Residual standard error: 25.18 on 6 degrees of freedom Multiple R-squared: 0.1109, Adjusted R-squared: -0.3336 F-statistic: 0.2495 on 3 and 6 DF, p-value: 0.8591
Given the Multiple R-squared value of this output, we can assume that this model is pretty awful. Regardless of such, our model could be effected by a phenomenon known as multiple collinearity. What this means, is that the variables could be interacting with each other, and thus, disrupting our model from providing accurate results. To prevent this from happening, we can measure the VIF() for each independent variable within the equation as it pertains to each other independent variable. If we were to do this manually, the code would resemble:
lm.multiregressx <- (lm(x ~ y + z )) # .4782 #
lm.multiregressy <- (lm(y ~ x + z )) # .5249 #
lm.multiregressz <- (lm(z ~ y + x )) # .1488 #
With each output, we would note the Multiple R-Squared value. I have provided those values to the right of each code line above. To then individually calculate the VIF(), we would utilize the following code for each R-Squared variable.
1 / (1 - .4282) # 1.916 #
1 / (1 - .5249) # 2.105 #
1 / (1 - .1488) # 1.175 #
This produces the values to the right of each code line. These values are the VIF() or Variance Influence Factor.
An easier way to derive the VIF() is to install the R package, “car”. Once “car” is installed, you can execute the following command:
vif(lm.multiregressw)
Which provides the output:
x y z 1.916437 2.104645 1.174747
Typically, most data scientists consider values of VIF() which are either over 5 or 10 (depending on sensitivity), to indicate that a variable ought to be removed. If you do plan on removing a variable from your model for such reasons, remove one variable at a time, as the removal of a single variable will effect subsequent VIF() measurements.
(Pearson) Coefficient of Correlation or cor():
Another tool at your disposal is the function cor(). Cor() allows the user to derive the coefficient of correlation ( r ) from numerical sets. For example:
As exciting as this seems, there is actually a better use for this function. For example, if you wanted to derive the coefficient of correlation as one variable as it pertains to others, you could use the following code lines to build a correlation matrix.
w x y z w 1.0000000 0.1057911 0.1836205 0.3261470 x 0.1057911 1.0000000 0.6914018 0.2546853 y 0.1836205 0.6914018 1.0000000 0.3853427 z 0.3261470 0.2546853 0.3853427 1.0000000
In the next article, we will again review the F-Statistic in preparation for a discussion pertaining to the concept of ANOVA. Until then, hang in there statistics fans!
Below is an example which illustrates a statistical concept:
# A new machine is purchased by a pastry shop which will be utilized to fill the insides of various pastries. A test run is commissioned in which the following proportions of jelly are inserted into each pastry. We are tasked with finding the 90% confidence interval for the mean filling amount from this process. #
s <- c(.28,.25,.22,.20,.33,.20)
w <- sd(s) / sqrt(length(s)) w
[1] 0.02092314
degrees <- length(s) - 1 degrees
[1] 5
t <- abs(qt(0.05, degrees)) t
[1] 2.015048
m <- (t * w) m
[1] 0.04216114
mean(s) + c(-m, m)
[1] 0.2045055 0.2888278
# We can state, with 90% confidence, that the machine will insert into pastries proportions of filling between .205 and .289. In the above execise, a new concept is introduced in step 4. The abs() function, in conjunction with the qt() function, allows you to find Critical-T values through the utilization of R without having to refer to a textbook. A few examples as to how this function can be utilized are listed below. #
# T VALUES (Confidence / Degrees of Freedom ) # abs(qt(0.25, 40)) # 75% confidence, 40 degrees of freedom, 1 sided (same as qt(0.75, 40)) # abs(qt(0.01, 40)) # 99% confidence, 40 degrees of freedom, 1 sided (same as qt(0.99, 40)) # abs(qt(0.01/2, 40)) # 99% confidence, 40 degrees of freedom, 2 sided #
I hope that you found this abbreviated article to be interesting and helpful. Until next time, I'll see you later Data Heads!
Welcome to the 60th article of Confessions of a Data Scientist. In this article we will be reviewing an incredibly important statistical topic: The T-Statistic.
The T-Statistic has many usages, but before we can begin to discuss these methods, I will briefly explain the conceptual significance of the subject matter, and when it is appropriate to utilize the T-Statistic.
The T-Statistic was created by a man named William Sealy Gosset. This concept was created by Gosset while he was working as a chemist at The Guinness Brewing Company. “Student”, was the anonymous name that was chosen by Gosset, which he utilized when publishing his findings. The reason for this anonymity, is that Guinness forbade its chemists from publishing their research. Gosset originally devised the T-Statistic and the concept of the “T-Test” as methods to measure the quality of stout. *
When to utilize The T-Statistic:
1. When the population standard deviation is unknown. (This is almost always the case).
2. When the current sample size is less than or equal to 30.
It is important that I mention now, that a common misconception is that the t-distribution assumes normality. This is not the case. Though, the t-distribution does begin to more closely resemble a normal distribution as the degrees of freedom approach 100. **
There are many uses for the t-distribution, some of which will be covered in future articles. However, the most commonly applied methods are: The One Sample T-Test, The Two Sample T-Test, and The Paired Test.
One Sample T-Test
This test is utilized to compare a sample mean to a specific value, it is used when the dependent variable is measured at the interval or ratio level.
To apply this test, you will need to define the following variables:
The size of the sample (N).
The applicable confidence interval, and the corresponding alpha.
The standard deviation of the sample (s).
The mean of the sample (M).
The mean of the population (mu).
To demonstrate the application of this method, I will provide two examples.
Example 1:
A factory employee believes that the cakes produced within his factory are being manufactured with excess amounts of corn syrup, thus altering the taste. 10 cakes were sampled from the most recent batch and tested for corn syrup composition. Typically, each cake should comprise of 20% corn syrup. Utilizing a 95 % confidence interval, can we assume that the new batch of cakes contain more than a 20% proportion of corn syrup?
The levels of the samples were:
.27, .31, .27, .34, .40, .29, .37, .14, .30, .20
Our hypothesis test will be:
H0: u = .2
HA: u > .2
With our hypothesis created, we can assess that mu = .2, and that this test will be right tailed.
We do do not to derive the additional variables, as R will perform that task for us. However, we must first input the sample data into R:
t.test(N, alternative = "greater", mu = .2, conf.level = 0.95)
# " alternative = " Specifies the type of test that R will perform. "greater" indicates a right tailed test. "left" indicates a left tailed test."two.sided" indicates a two tailed test. #
Which produces the output:
One Sample t-test
data: N t = 3.6713, df = 9, p-value = 0.002572 alternative hypothesis: true mean is greater than 0.2 95 percent confidence interval: 0.244562 Inf sample estimates: mean of x 0.289
t = The t-test statistic, which derived by hand would be: the sample mean (M), minus the population mean (mu), divided by the standard deviation (S) divided by the square root of the size of the sample (N).
df = This value is utilized in conjunction with the confidence interval to define critical-t. It is derived by subtracting 1 from the value of N.
p-value = This is the probability that, given the parameters of the test, that we will obtain a value equal or greater than the t-region defined within our model.
##% confidence interval = This describes the numerical values which exists on each side of the t-distribution, which between both, contains ##% of the values of the model.
sample estimates = This value is the mean of the sample(s) utilized within the function.
With this output we can conclude:
With a p-value of .002572 (.002572 < .05), and a corresponding t-value of 3.6713, we can conclude that, at a 95% confidence interval, that the cakes are being produce contain an excess amount of corn syrup.
Example 2:
A new type of gasoline has been synthesized by a company that specializes in green energy solutions. The same fuel type, when synthesized traditionally, typically provides coup class vehicles with a mean of 25.6 mpg. Company statisticians were provided with a test sample of 8 cars which utilized the new fuel type. This sample of vehicles utilized the fuel at a mean rate of 23.2 mpg, with a sample standard deviation of 5 mpg. Can it be assumed, at a 95% confidence interval, that vehicle performance significantly differs when the new fuel type is utilized?
Our hypothesis will be:
H0: u = 25.6
HA: u NE 25.6
With our hypothesis created, we can assess that mu = 25.6, and that this test will be two tailed.
In this case we will need to create a random sample within R that adheres to the following parameters.
N = 8 Sample Size
S = 5 Sample Standard Deviation
M = 23.2 Sample Mean
This can be achieved with the following code block: sampsize <- 8 standarddev <- 5 meanval <- 23.2
N <- rnorm(sampsize) N <- standarddev*(N-mean(N))/sd(N)+meanval With our “N” variable defined, we can now let R do the rest of the work:
t.test(N, alternative = "two.sided", mu = 25.6, conf.level = 0.95)
Which produces the output:
One Sample t-test
data: N t = -1.3576, df = 7, p-value = 0.2167 alternative hypothesis: true mean is not equal to 25.6 95 percent confidence interval: 19.0199 27.3801 sample estimates: mean of x 23.2
With this output we can conclude:
With a p-value of 0.2167 (.2167 > .05), and a corresponding t-value of -1.3576, we cannot conclude that, at a 95% confidence interval, that vehicle performance differs depending on utilization of fuel type.
Two Sample T-Test
This test is utilized if you randomly sample different sets of items from two separate control groups.
Example:
A scientist creates a chemical which he believes changes the temperature of water. He applies this chemical to water and takes the following measurements:
70, 74, 76, 72, 75, 74, 71, 71
He then measures temperature in samples which the chemical was not applied.
74, 75, 73, 76, 74, 77, 78, 75
Can the scientist conclude, with a 95% confidence interval, that his chemical is in some way altering the temperature of the water?
For this, we will use the code: N1 <- c(70, 74, 76, 72, 75, 74, 71, 71)
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.3007929 4.4492071
sample estimates:
mean of x mean of y
75.250 72.875
# Note: In this case, the 95 percent confidence interval is measuring the difference of the mean values of the samples. #
# An additional option is available when running a two sample t-test, The Welch Two Sample T-Test. To utilize this option while performing a t-test, the "var.equal = TRUE" must be changed to "var.equal = FALSE". The output produced from a Welch Two Sample t-test is slightly more robust and accounts for differing sample sizes. #
From this output we can conclude:
With a p-value of 0.02773 (.0.02773 < .05), and a corresponding t-value of 2.4558, we can state that, at a 95% confidence interval, that the scientist's chemical is altering the temperature of the water.
Paired T-Test
This test is utilized if you are sampling the same set twice, once for each variable.
Example:
A watch manufacturer believes that by changing to a new battery supplier, that the watches that are shipped which include an initial battery, will maintain longer lifespan. To test this theory, twelve watches are tested for duration of lifespan with the original battery.
The same twelve watches are then re-rested for duration with the new battery.
Can the watch manufacturer conclude, that the new battery increases the duration of lifespan for the manufactured watches? (We will assume an alpha value of .05).
For this, we will utilize the code: N1 <- c(376, 293, 210, 264, 297, 380, 398, 303, 324, 368, 382, 309) N2 <- c(337, 341, 316, 351, 371, 440, 312, 416, 445, 354, 444, 326)
t.test(N2, N1, alternative = "greater", paired=TRUE, conf.level = 0.95 )
This produces the output: Paired t-test
data: N2 and N1 t = 2.4581, df = 11, p-value = 0.01589 alternative hypothesis: true difference in means is greater than 0 95 percent confidence interval: 12.32551 Inf sample estimates: mean of the differences 45.75
From this output we can state:
With a p-value of 0.01589 (0.01589 < .05), and a corresponding t-value of 2.4581, we can conclude that, at a 95% confidence interval, that the new battery increases the duration of lifespan for the manufactured watches.
Below are two examples which illustrate a statistical concept:
The Distribution of Sample Means:
The average mortgage for new families in Texas in $ 1300, with a standard deviation of $ 800. If 100 new families are surveyed, what is the probability that the mean mortgage payment will exceed $ 1500?
800/sqrt(100) [1] 80
(1500-1300) / 80
[1] 2.5
pnorm(2.5, lower.tail=FALSE)
[1] 0.006209665
There is a 0.62% chance of this occurring.
An American football team typically scores 17 points per game, with a standard deviation of 3 points. In a sample of 16 games, what is the probability that the team scored an average of between 16 - 19 points? 3/sqrt(16)
