DESCRIPTIVE ANALYSIS AND VISUALIZATION IN R USING PRESTIGE DATASET IN THE CAR PACKAGE

Below is a warm-up graph to get you going with R.

plot(cars, col = "darkblue", cex = 2)

The dataset used for this analysis is Prestige from the package car. First, we need to install the required packages for our analysis.

# Load the packages that contain the dataset and the data viz package
library(car) # Package containing the dataset
library(ggplot2) # Package for some nice and cool visualizations
library(MASS) # Package for box-cox transformation

Prestige is a dataset of Canadian Occupations. The Prestige data frame contains 102 rows and 6 columns. The observations relate to occupations. The data frame includes the following columns:

1. Variable education is the average education of occupational incumbents, years, in 1971.

2. variable income is the average income of incumbents, in 1971.

3. variable women is the percentage of incumbents who are women.

4. Variable prestige is Pineo-Porter score for occupation, from a social survey conducted in mid-1960s.

5. Variable census is the Canadian Census occupational code.

6. Variable type relates to the type of occupation. A factor with levels (note: out of order): bc, Blue Collar; prof, Professional, Management, and Technical; wc, White Collar.

Source of the dataset is: Canada(1971) Census of Canada. Vol. 3, Part 6. Statistics Canada [pp. 19-1-19-21].

1. SUMMARY and DESCRIPTIVE STATISTICS

Summary (or descriptive) statistics are the first figures used to represent nearly every dataset. They also form the foundations for more complicated computations and analyses. Therefore, they are essential to the analysis process. In this paper, we will explore the ways in which R can be used to calculate summary statistics, including mean, standard deviation, range, and percentile. Also included here is the summary function, which is one of the most useful tools in the R set of commands.

First, let us inspect the Prestige dataset.

# Display or print the first 10 observations of our dataset
print(head(Prestige, n = 10))

We observe that the rows of our dataset refer to occupations. Some of the occupations are gov.administrators, general.managers, accountants, chemists, physicists, biologists, architects, civil.engineers and mining.engineers etc. In addition, each occupation has its own records. We can also explore the nature of variables/columns we have in the Prestige dataset, using ls(DATAVAR) Or names(DATAVAR).

# To see the variable in the dataset; Use names(dataset) or ls(dataset)
ls(Prestige)

[1] "census"    "education" "income"    "prestige"  "type"      "women"    

# To see the number of columns and number of rows in the Prestige dataset; use ncol(dataset) and nrow(dataset)
ncol(Prestige)

[1] 6

nrow(Prestige)

[1] 102

From the result output, the Prestige dataset contains 6 variables and 102 rows. A more advanced technique to see the structure of the dataset is to use the str(DATAVAR) function.

# A more advanced and complete way to see the structure of our dataset
str(Prestige)

'data.frame':   102 obs. of  6 variables:
 $ education: num  13.1 12.3 12.8 11.4 14.6 ...
 $ income   : int  12351 25879 9271 8865 8403 11030 8258 14163 11377 11023 ...
 $ women    : num  11.16 4.02 15.7 9.11 11.68 ...
 $ prestige : num  68.8 69.1 63.4 56.8 73.5 77.6 72.6 78.1 73.1 68.8 ...
 $ census   : int  1113 1130 1171 1175 2111 2113 2133 2141 2143 2153 ...
 $ type     : Factor w/ 3 levels "bc","prof","wc": 2 2 2 2 2 2 2 2 2 2 ...

Prestige is a data.frame, a datatype with more than one row and column. Prestige includes 5 numeric variables and one categorical variable. From here, we can now perform a summary and descriptive statistics of our dataset.

MEAN of EACH VARIABLE

In R, a mean can be calculated on an isolated variable via the mean(VAR) command, VAR is the name of the variable whose mean we want to compute. Alternative, the mean can be calculated for all the variables in the dataset using mean(DATAVAR) function, where DATAVAR is the name of the dataset containing the variables. For analysis purposes, we are going to exclude the variables census and type from the descriptive statistics. To select a subset of a dataset, use subset(DATAVAR, select = c(“VAR1”, “VAR2”, “VAR3”….“VARi”)) command. You can also type ?subset() in your R console followed by ENTER to learn how to subset() vectors, matrices and data frames. The code below demonstrates how to select a subset of Prestige dataset.

# Subsetting education, income, women, and prestige from the dataset Prestige
subset.data <- subset(Prestige, select = c("education", "income", "women", "prestige"))

# Checking subset.data to make sure we have the needed subset
str(subset.data)

'data.frame':   102 obs. of  4 variables:
 $ education: num  13.1 12.3 12.8 11.4 14.6 ...
 $ income   : int  12351 25879 9271 8865 8403 11030 8258 14163 11377 11023 ...
 $ women    : num  11.16 4.02 15.7 9.11 11.68 ...
 $ prestige : num  68.8 69.1 63.4 56.8 73.5 77.6 72.6 78.1 73.1 68.8 ...

From the above output, we see that we got the subset we need. Let’s find the mean of each variable in the selected subset.

# Calculate the mean of a variable with mean(DATAVAR$VAR); mean of variable education
mean(subset.data$education)

[1] 10.73804

The mean year of education was 11 years in Canada in 1971. That means that the average Canadian had a least a High School education in 1971.

# Mean of income variable
mean(subset.data$income)

[1] 6797.902

Here the mean income for incumbents in 1971 is about $6,798.

# Mean of percentage of incumbents who are women
mean(subset.data$women)

[1] 28.97902

The average percentage of incumbents in the survey who are women is about 29%.

# The last variable is prestige
mean(subset.data$prestige)

[1] 46.83333

The mean score of Pineo_Porter measure for occupation is about 47. The computations above are done on individual variables. Alternatively, a mean can be calculated for each of the variables in the dataset by using the formula summary(dataset). We will approach this alternative later. Let’s find the standard deviation of each variable in our subset dataset.

STANDARD DEVIATION OF EACH VARIABLE

Standard deviations are calculated in the same way as means within R. The standard deviation of a single variable can be computed using the formula sd(VAR), where VAR is the name of the variable whose standard deviation you want to retrieve. Standard deviation measures how spread your data are. The codes below demonstrate the use of the standard deviation function.

# What is the standard deviation of education?
sd(subset.data$education)

[1] 2.728444

# What is the standard deviation of income?
sd(subset.data$income)

[1] 4245.922

# What is the standard deviation of women?
sd(subset.data$women)

[1] 31.72493

# What is the standard deviation of prestige?
sd(subset.data$prestige)

[1] 17.20449

RANGE of EACH VARIABLE : MINIMUM & MAXIMUM

Continuing in the same trajectory, minimum can be computed on a single variable using the min(VAR) formula. In the same token, max(VAR) operates similarly. Minimum and Maximum give the min and max of individual variables in the dataset. The codes below show how to calculate minimums and maximums.

# Minimum and maximum years of education of survey participants in 1971
min(subset.data$education) 

[1] 6.38

max(subset.data$education)

[1] 15.97

The minimum year of incumbents’ education in 1971 was about 6 years. And the maximum year of participants’ education was 16 years, meaning we had some college-educated incumbents in the survey.

# Minimum and Maximum years of income of survey participants in 1971
min(subset.data$income); max(subset.data$income)

[1] 611
[1] 25879

From the output, the minimum income is 611 and the maximum income is about $26000. This shows large gaps in income distribution among survey participants. The maximum income is about four times the minimum one.

# Minimum and Maximum percentage of incumbents who are women.
min(subset.data$women); max(subset.data$women)

[1] 0
[1] 97.51

From the output, we may infer that a large percentage of women participated in the survey.

# Minimum and maximum score for occupation
min(subset.data$prestige) ; max(subset.data$prestige)

[1] 14.8
[1] 87.2

The range for prestige scores is from ~ 15 to ~ 87. The Range() function can provide minimums and maximums of different variables in our subset dataset.

RANGE

The range of a particular variable, that is, its minimum and maximum, can be retrieved using the range(VAR) command. Like with min and max functions, using range(dataset) is not very useful since it considers the entire dataset, rather than each individual variable. Consequently, it is recommended that ranges be computed on individual variables. These computations are demonstrated in the following codes:

# Calculate the range of a variable with range(VAR)
# Range of variable education
range(subset.data$education)

[1]  6.38 15.97

# Range of variable income
range(subset.data$income)

[1]   611 25879

# Range of variable women
range(subset.data$women)

[1]  0.00 97.51

# Range of variable prestige
range(subset.data$prestige)

[1] 14.8 87.2

PERCENTILES : VALUES from PERCENTILES (QUANTILE)

You can get more insight into the distribution of a set of observations by examining quantiles. A quantile is a value computed from a collection of numeric measurements that indicates an observation’s rank when compared to all other present observations. Alternatively, quantile can be expressed as a percentile, this is identical but on a percent scale of 0 to 100.

Obtaining quantiles and percentile in R is done using the quantile() function. The command is written quantile(VAR, c(PROB1, PROB2, PROB3,….PROBi)) or quantile(VAR, prob = c(prob value1, prob value 2, prob value 3…prob valuei)).

# Calculate the 25th, 50th, 75th, and 95th percentiles for education in the subset dataset
quantile(subset.data$education, prob = c(.25, .50, .75, .95))

    25%     50%     75%     95% 
 8.4450 10.5400 12.6475 15.4290 

quantile(subset.data$education, c(.25, .50, .75, .95))

    25%     50%     75%     95% 
 8.4450 10.5400 12.6475 15.4290 

We have used two slightly different codes that lead to the same output. From the output, we see that 25% of average years of education was about 8 years of education, the median year of education was about 10.5 years of education, 75% of average education of incumbents was about 13 years of education, and 95% of participants have about 16 years of education. Also, this may explain the wide gap in wage range since education seems to be correlated to income.

# Calculate the 25th, 50th, 75th, and 95th percentiles for income in the dataset
quantile(subset.data$income, prob = c(.25, .50, .75, .95))

     25%      50%      75%      95% 
 4106.00  5930.50  8187.25 14156.45 

For the income variable, let find the 0 th, the 25th , the 50th, the 75th, and the 100th percentiles.

# Calculate the 0th, 25th, 50th, 75th, and 100th percentiles
quantile(subset.data$income, c(0, .25, .50, .75, 1))

      0%      25%      50%      75%     100% 
  611.00  4106.00  5930.50  8187.25 25879.00 

The median income of survey participants was $US 5930.50 in 1971. Higher income tended to be in the upper bounds of the distribution. Income distribution seems to reflect the levels of education of participants in the survey. The more educated you are, the higher the likelihood of being in the upper income’s bound.

# Calculate the 0th, 25th, 50th, 75th, and 100th perentiles for women
quantile(subset.data$women, prob = c(0, .25, .50, .75, 1))

     0%     25%     50%     75%    100% 
 0.0000  3.5925 13.6000 52.2025 97.5100 

The minimum proportion of incumbents who are women is 0, and the maximum proportion of participants who are women is ~ 98. Finally, the median number of women is ~ 17.

# Calculate the 0th, 25th, 50th, 75th, and 100th percentiles for prestige
quantile(subset.data$prestige, c(0, .25, .50, .75, 1))

    0%    25%    50%    75%   100% 
14.800 35.225 43.600 59.275 87.200 

Here we have used the quantile function to obtain what is called the FIVE-NUMBER SUMMARY of subset.data, comprised of 0th percentile(the minimum), the 25th percentile, the 50th percentile, and the 100th percentile(the maximum). The 0.25th quantile is referred to as the first or LOWER QUARTILE, and 0.75th quantile is referred to as the third or UPPER QUARTILE. Also, the 50th percentile refers to the MEDIAN. The median is the second quartile, with the maximum value being the fourth quartile.

PERCENTILES FROM VALUES (PERCENTILE RANK)

In the opposite situation, where a percentile rank corresponding to a given value is needed, one has to devise a custom method. Here are the steps involved in computing a percentile rank.

1. Count the number of data points that are at or below the given value
2. Divide by the total number of data points
3. multiple by 100

The formula for calculating a percentile rank can be derived from this command: percentile rank = length(VAR[VAR <= VALUE]) / length(VAR) * 100. length(VAR[VAR <= VALUE]) counts the number of data points in a variable that are below the given value. The ‘<=’ can be replaced with other operators, such as ‘<’, ‘>’, and =. The length(VAR) counts the number of data points in the variable. The final step is to multiply the result by 100 to transform the decimal value into a percentage. Let’s apply these steps in the following examples:

# In the sample, 8 years of education is at what percentile rank?
length(subset.data$education[subset.data$education <= 8]) / length(subset.data$education) * 100

[1] 19.60784

8 years of education is at 19.60 percentile rank.

# In the sample, $10000 is at what percentile rank?
length(subset.data$income[subset.data$income <= 10000]) / length(subset.data$income) * 100

[1] 87.2549

# In the sample, $10000 is at what percentile rank?
length(subset.data$income[subset.data$income >= 10000]) / length(subset.data$income) * 100

[1] 12.7451

The percentile rank reveals some interesting facts on income distribution of survey incumbents. About 13 % of participants earned more than or equal to $10,000.00 in 1971, while about 87 % earned less than or equal to 10,000.00. Similarly, let’s find the percentile rank of income above and below, say 15,000 dollars.

# In the sample, income less than or equal to say 15000 is at what percentile rank?
length(subset.data$income[subset.data$income <= 15000]) / length(subset.data$income) * 100

[1] 96.07843

# In the sample, income greater than or equal to say 15000 is at percentile rank?
length(subset.data$income[subset.data$income >= 15000]) / length(subset.data$income) * 100

[1] 3.921569

These two outputs confirm the wide gap in income distribution among survey participants. A useful multipurpose function in R is the use of the SUMMARY(X) function, as I previously mentioned.

SUMMARY

A very useful function in R is summary(x), where x can be one of any number of objects, including datasets, variables, and linear models. When used, the summary(x) provides summary data related to the individual object that is included into it. Thus, the summary function has a different output depending on what kind of object it takes as an argument. This method is valuable because it often sums up what we previously did and provides exactly what is needed in summary statistics. Let’s apply this command to the sample dataset.

# Summarize a variable using summary(VAR). Summary statistics of education
print(summary(subset.data$education))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  6.380   8.445  10.540  10.740  12.650  15.970 

We get the min, 1st quartile, Median, Mean, 3rd quartile, and the maximum years of education. This step provides more summary statistics information. We will apply,the same command this time,to the dataset.

# Summarize the subset.data sample using the command summary(subset.data)
print(summary(subset.data))

   education          income          women           prestige    
 Min.   : 6.380   Min.   :  611   Min.   : 0.000   Min.   :14.80  
 1st Qu.: 8.445   1st Qu.: 4106   1st Qu.: 3.592   1st Qu.:35.23  
 Median :10.540   Median : 5930   Median :13.600   Median :43.60  
 Mean   :10.738   Mean   : 6798   Mean   :28.979   Mean   :46.83  
 3rd Qu.:12.648   3rd Qu.: 8187   3rd Qu.:52.203   3rd Qu.:59.27  
 Max.   :15.970   Max.   :25879   Max.   :97.510   Max.   :87.20  

The output of the preceding summary provides the descriptive statistics of all objects in the sample data set. Under each variable, we have its summary statistics. Now that we know how to do summary statistics, we can delve into the visual part of the analysis to see the relation between the variables. The visual part will be done using basic R graphics packages and ggplot2, thanks to Hadley Wickham, Chief Scientist at Rstudio.

DATA VISUALIZATION

1. MATRIX of PLOTS

The single type of planar scatterplot is really useful only when comparing two numeric-continuous variables. When there are more continuous variables of interest, it is possible to display this information satisfactorily on a single plot. A simple and common solution is to generate a two-variable scatterplot for each pair of variables and show them together in a structured way; this is referred to as a Scatterplot Matrix. We have four continuous/numeric variables in the subset dataset we have selected. Working with base R graphics, use the pairs function.

# Drawing a scatterplot matrix of education, income, women, and prestige using the pairs function
pairs(subset.data, pch = 16, col = "blue", main = "Matrix Scatterplot of Education, Income, Women, Prestige")

Our matrix scatterplot may be too big to fit on our screen. However, if you run the above code in your console, you will get the matrix plot in Rstudio’s graphics area. The interpretation of the above plots depends upon the labeling of the diagonal panels, running from the top left to the bottom right. They will appear in the same order as the columns given as the first argument. These “label panels” allow you to determine which individual plot in the matrix corresponds to each pair of variables. For instance, the first column of the scatterplot matrix corresponds to an x-axis variable of education; the third row of the matrix corresponds to a y-axis variable of women, and each row and column displays a scale that is constant moving left/right or up/down, respectively. The plot of prestige(y) on income(x) at position row 4 and column 2 displays the same data as the scatterplot at position row 2 and column 4 but flipped on its axis. Likewise, the plot of income(y) on education(x) at position row 2 and column 1 displays the same data as the scatterplot at column 2 row 1 flippedd on its axis. The scatterplot matrices therefore allow for easier comparison of pairwise relationships formed by observations made on multiple continuous variables. Instead of viewing a panoramic relationship between all the variables in our subset dataset, let’s use simple scatter plots to visualize the relationship between two variables.

2. SCATTERPLOT

A scatterplot is a useful way to visualize the relationship between two variables displayed as x-y coordinate plots. Similar to correlations, scatterplots are often used to make initial diagnoses before any statistical analyses are performed.

The simplest way to create a scatterplot is to directly graph two variables using the default settings. In R, this can be achieved using the formula plot(VARX, VARY) function, where VARX is the variable to plot along the x-axis and VARY is the variable to plot along the y-axis. I will also add the ggplot2 version for the scatterplot. Let’s look at the relationship between education and income.

# Picture of the relationship between education and income
plot(subset.data$education, subset.data$income, xlab = "Education", ylab = "Income", pch = 16, col = "blue",
     main = " Scatterplot of Education vs. Income") 

# Same plot nicely done with ggplot2
ggplot(subset.data) +
  geom_point(aes(x = education, y = income), col = 'blue', size = 3) +
  ggtitle("Education vs. Income Scatterplot") +
  theme(plot.title = element_text(hjust = 0.5))

From the output, we see that as the average of years of education and income are positively related. That is higher education is rewarded by higher income. The scatterplot also displays the presence of potential outliers. We can further expand our visual analysis by calculating the slope and intercept of line of best fit.

# Calculate slope and intercept of line of best fit
coef(lm(income ~ education, data = subset.data))

(Intercept)   education 
 -2853.5856    898.8128 

We will use geom_abline() function to add the intercept avlue and estimated coefficient of education to our plot.

# Adding a line of best fit; intercept and slope
ggplot(subset.data) +
  
  geom_point(aes(x = education, y = income), col = "blue", size = 3) +
  
  geom_abline(aes(intercept = -2853, slope = 899), col = "darkred") +
  
  ggtitle("Education vs. Income Scatterplot With The Best Fit Line") +
  
  theme(plot.title = element_text(hjust = 0.5))

We can use our original Prestige dataset to display the income by types of occupation. In our dataset, we have type, a factor variable, that refers to three occupational levels: bc for Blue Collar; prof for Professional, Managerial, and Technical; and wc for White Collar. Let’s construct a bar plot for this categorial variable.

# Easy way to do it using the plot() function
plot(Prestige$type, main = "Bar Plot of Occupational Types", col = "cyan", ylab = "COUNT")

The output demonstrates that there are more blue collar workers in the survey than any other occupational groups. That may explain the disproportional distribution of income.

# Same plot using ggplot2
ggplot(Prestige) +
  
  geom_bar(aes(x = type)) +
  
  ggtitle("Count of Incumbents' Occupation") +
  
  theme(plot.title = element_text(hjust = 0.5))

We clearly see the difference between the two plots. The first one doesn’t take into account NAs occupation type, but the second plot using ggplot takes into account the NAs. Moving forward, we are going to plot income against types of occupations and add labels to these occupations. We are going to install a useful package called ggrepel that labels the scatterplot’s points.

ggplot(Prestige, aes(x = education, y = income)) +
  
  geom_text(aes(label = type), size = 4) +
  
   ggtitle("Scatterplot With Data Labels") +
  
  theme(plot.title = element_text(hjust = 0.5))

From the output, we see that Professionals have higher average years of education and earn higher income than other groups. Blue collars have less average years of education and are in the lower bounds of income.

# Using ggrepel to come up with the same plot
install.packages("ggrepel")

Installing package into <U+393C><U+3E31>C:/Users/pkolo/Documents/R/win-library/3.3<U+393C><U+3E32>
(as <U+393C><U+3E31>lib<U+393C><U+3E32> is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.3/ggrepel_0.6.5.zip'
Content type 'application/zip' length 936513 bytes (914 KB)
downloaded 914 KB

package ‘ggrepel’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\pkolo\AppData\Local\Temp\RtmpCw3QNq\downloaded_packages

library(ggrepel)
ggplot(Prestige, mapping = aes(x = education, y = income)) +
  
  geom_point() +
  
  geom_text_repel(aes(label = type), size = 4) +
  
  ggtitle("Scatterpot With Data Points and Labels") +
  
  theme(plot.title = element_text(hjust = 0.5))

We can map variables to other variables as well.

ggplot(Prestige, mapping = aes(x = education, y = income)) +
  
  geom_point(aes(color = prestige, shape = type)) +
  
  ggtitle("Education vs. Income Scatterpot Per Prestige and Type") +
  
  theme(plot.title = element_text(hjust = 0.5))

In this output we map income and education with type and prestige. The output demonstrates that blue collars and white collar workers have low prestige scores, compared to professional. Income and levels of education seem to explain these trends.