In this report, the Prostate cancer dataset will be used for analysis. Comparison of the data will be analyzed by using the Tree and Regression analyses for predicting both lpsa and lcavol. In total, there will be four separate analyses that you compare. The prostate cancer dataset is a study by Stamey et al. (1989) that examined the correlation between the level of prostate specific antigen (PSA) and a number of clinical measures, in 97 men who were about to receive a radical prostatectomy.
Install and Load Packages
There will be multiple libraries we need to install and load in this R tutorial before we load the data.
Input:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | install.packages('C50') install.packages('gmodels') install.packages('party') install.packages('RColorBrewer') install.packages('RWeka') install.packages('psych') install.packages('rpart') install.packages('rpart.plot') library(C50) library(gmodels) library(party) library(RColorBrewer) library(RWeka) library(psych) library(rpart) library(rpart.plot) |
Since we will be using the concrete data set, you will need to download this data set. This dataset is already packaged and available for an easy download from the dataset page. or directly from here Prostate – prostate.csv
Load and View the Prostate Data
In the below, we will load the prostate data file and view the contents of the data with the summary() function and build histrograms with the histogram() function.
summary() function
Input:
1 2 3 | prostate <- read.csv("prostate.csv", stringsAsFactors = TRUE) summary(prostate) |
Output:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | lcavol age lbph lcp gleason Min. :-1.3471 Min. :41.00 Min. :-1.3863 Min. :-1.3863 Min. :6.000 1st Qu.: 0.5128 1st Qu.:60.00 1st Qu.:-1.3863 1st Qu.:-1.3863 1st Qu.:6.000 Median : 1.4469 Median :65.00 Median : 0.3001 Median :-0.7985 Median :7.000 Mean : 1.3500 Mean :63.87 Mean : 0.1004 Mean :-0.1794 Mean :6.753 3rd Qu.: 2.1270 3rd Qu.:68.00 3rd Qu.: 1.5581 3rd Qu.: 1.1787 3rd Qu.:7.000 Max. : 3.8210 Max. :79.00 Max. : 2.3263 Max. : 2.9042 Max. :9.000 lpsa Min. :-0.4308 1st Qu.: 1.7317 Median : 2.5915 Mean : 2.4784 3rd Qu.: 3.0564 Max. : 5.5829 |
hist() function
Input:
1 | hist(prostate$lpsa) |
Output:
Input:
1 | hist(prostate$lcavol) |
Output:
Rplot
This classification is comparing the variable of type to all predictors within the prostate dataset. We will be comparing two predictors; lcavol and lpsa.
LCAVOL Rplot
Input:
1 2 3 | lcavol_modelR <- lm(lcavol ~ age + lbph + lcp + lpsa + gleason, data = prostate) lcavol_modelR |
Output:
1 2 3 4 5 6 | Call: lm(formula = lcavol ~ age + lbph + lcp + lpsa + gleason, data = prostate) Coefficients: (Intercept) age lbph lcp lpsa gleason -1.49371 0.01902 -0.08918 0.29727 0.53955 0.05240 |
LPSA Rplot
Input:
1 2 3 | lpsa_modelR <- lm(lpsa ~ lcavol + age + lbph + lcp + gleason, data = prostate) lpsa_modelR |
Output:
1 2 3 4 5 6 | Call: lm(formula = lpsa ~ lcavol + age + lbph + lcp + gleason, data = prostate) Coefficients: (Intercept) lcavol age lbph lcp gleason 1.83165 0.65523 -0.01027 0.14506 0.07173 0.06168 |
Train and Test Data
rpart Formula Implementation
We can specify age as the outcome variable and use the dot notation to allow all the other columns in the credit_train data frame to to be used as predictors.
Input:
1 2 3 4 5 | lcavol_pred <- predict(lcavol_modelR, prostate) lpsa_pred <- predict(lpsa_modelR, prostate) prostate_train <- prostate[1:60, ] prostate_test <- prostate[61:97, ] |
LCAVOL Train Data
Input:
1 2 3 | lcavol.rpart <- rpart(lcavol ~ age + lbph + lcp + gleason, data = prostate_train) lcavol.rpart |
Output:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | n= 60 node), split, n, deviance, yval * denotes terminal node 1) root 60 61.667440 0.81724440 2) lcp< 0.2101769 50 41.521730 0.57621020 4) gleason< 6.5 28 22.037790 0.27570190 8) age< 63.5 17 12.795910 -0.08539058 * 9) age>=63.5 11 3.599646 0.83375390 * 5) gleason>=6.5 22 13.737240 0.95867530 10) age>=63.5 15 7.008366 0.73164000 * 11) age< 63.5 7 4.298893 1.44518000 * 3) lcp>=0.2101769 10 2.716462 2.02241500 * |
LPSA Train Data
Input:
1 2 3 | lpsa.rpart <- rpart(lpsa ~ age + lbph + lcp + gleason, data = prostate_train) lpsa.rpart |
Output:
1 2 3 4 5 6 7 8 9 10 11 12 | n= 60 node), split, n, deviance, yval * denotes terminal node 1) root 60 38.815370 1.809756 2) gleason< 6.5 30 20.284900 1.481545 4) lbph< 0.9410666 19 12.333000 1.192772 * 5) lbph>=0.9410666 11 3.630779 1.980336 * 3) gleason>=6.5 30 12.067120 2.137967 6) lbph< 0.3306992 14 7.923121 1.927040 * 7) lbph>=0.3306992 16 2.976120 2.322529 * |
Decision Trees
As one can see below, these decision trees are a big improvement from the earlier chapter. This gives greater detail on the outcome and predictors for the credit set. I prefer the rplot with the color coded decision trees for percentages and yes and no outcomes. Let’s visualize the data for lcavol and lpsa.
LCAVOL Decision Tree
Input:
1 | rpart.plot(lcavol.rpart, digits = 4, fallen.leaves = TRUE, type = 3, extra = 101) |
Output:
LPSA Decision Tree
Input:
1 | rpart.plot(lpsa.rpart, digits = 4, fallen.leaves = TRUE, type = 3, extra = 101) |
Output:
Model Performance Evaluation
LCAVOL
Input:
1 2 3 | p.lcavol <- predict(lcavol.rpart, prostate_test) summary(p.lcavol) |
Output:
1 2 | Min. 1st Qu. Median Mean 3rd Qu. Max. -0.08539 1.44518 2.02242 1.66620 2.02242 2.02242 |
LPSA
Input:
1 2 3 | p.lpsa <- predict(lpsa.rpart, prostate_test) summary(p.lpsa) |
Output:
1 2 | Min. 1st Qu. Median Mean 3rd Qu. Max. 1.193 1.927 1.980 2.063 2.323 2.323 |
MAE
Measuring correlation on lcavol and lpsa.
LCAVOL
Input:
1 | cor(p.lcavol, prostate_test$lcavol) |
Output:
1 | [1] 0.3837933 |
LPSA
Input:
1 | cor(p.lpsa, prostate_test$lpsa) |
Output:
1 | [1] 0.05227984 |
Actual and Predicted function for MAE
1 2 3 | MAE <- function(actual, predicted) { mean(abs(actual - predicted)) } |
LCAVOL MAE
Input:
1 2 3 4 5 | MAE(p.lcavol, prostate_test$lcavol) mean(prostate_test$lcavol) MAE(2.21, prostate_test$lcavol) |
Output:
1 2 3 4 5 | [1] 0.7968151 [1] 2.213953 [1] 0.6830092 |
LPSA MAE
Input:
1 2 3 4 5 | MAE(p.lpsa, prostate_test$lpsa) mean(prostate_test$lpsa) MAE(3.56, prostate_test$lpsa) |
Output:
1 2 3 4 5 | [1] 1.499867 [1] 3.562653 [1] 0.5247962 |
Regression
Correlation matrix for age, lbph, lcp, gleason
Input:
1 | cor(prostate[c("lcavol", "age", "lbph", "lcp", "gleason", "lpsa")]) |
Output:
1 2 3 4 5 6 7 | lcavol age lbph lcp gleason lpsa lcavol 1.0000000 0.2249999 0.027349703 0.675310484 0.43241706 0.7344603 age 0.2249999 1.0000000 0.350185896 0.127667752 0.26889160 0.1695928 lbph 0.0273497 0.3501859 1.000000000 -0.006999431 0.07782045 0.1798094 lcp 0.6753105 0.1276678 -0.006999431 1.000000000 0.51483006 0.5488132 gleason 0.4324171 0.2688916 0.077820447 0.514830063 1.00000000 0.3689868 lpsa 0.7344603 0.1695928 0.179809410 0.548813169 0.36898680 1.0000000 |
scatterplot matrix using pairs()
Input:
1 | pairs(prostate[c("lcavol", "age", "lbph", "lcp", "gleason", "lpsa")]) |
Output:
Correlation Ellipse
The oval-shaped object on each scatterplot is called a correlation ellipse. This shows a visual of the correlation strength. The dot at the center of the correlation ellipse, indicates the point at the mean values for the x and y axis variables. The correlation between the two variables is indicated by the shape of the ellipse is shown by the shape more stretched, the stronger the correlation. When there’s a perfect rounded oval, this indicated a very weak correlation.
Input:
1 | pairs.panels(prostate[c("lcavol", "age", "lbph", "lcp", "gleason", "lpsa")]) |
Output:
Implementing lm() function
The lm() function will return a regression model object that can be used to make predictions. Interactions between independent variables can be specified using the * operator.
Dependent lcavol against all other variables
1 2 3 | > pro_model1 <- lm(lcavol ~ age + lbph + lcp + gleason, data = prostate) > pro_model1 |
1 2 3 4 5 6 | Call: lm(formula = lcavol ~ age + lbph + lcp + gleason, data = prostate) Coefficients: (Intercept) age lbph lcp gleason -0.78185 0.02085 -0.01689 0.51970 0.13253 |
Dependent lpsa against all other variables
1 2 3 | > pro_model2 <- lm(lpsa ~ age + lbph + lcp + gleason, data = prostate) > pro_model2 |
1 2 3 4 5 6 | Call: lm(formula = lpsa ~ age + lbph + lcp + gleason, data = prostate) Coefficients: (Intercept) age lbph lcp gleason 1.319355 0.003393 0.133996 0.412253 0.148513 |
Improving Model Performance
The below break down for the lm() and using summary() is showing residuals and coefficients. The residuals section provides a statistical summary for the prediction errors. Some of these errors are quite substantial. As one can see, 50 percent of the errors fall into the 1Q and 3Q values. For the dependent variable lcavol, this means that predictions were between -0.527 over the true value and 0.578 were under the true value. While lpsa were between -0.444 over the true value and 0.522 were under the true value
For each of the estimated regression coefficient, the p-value provides an estimate of the probability that the true coefficient is given zero given the value of the estimate.The below model we created, has several highly significant variables, and they seem to be related to the outcomes in logical ways.
The coefficient of determination, also known as the multiple R-squared value, provides a measure of how well our model as a whole explains the values of the dependent variable. This means that the closer the value is to 1.0, the better the model explains the data. As one can see the below data, the adjusted R-squared data is 0.664 for lcavol and 0.575 for lpsa. This shows us that the 66 percent of the variation is in the dependent variable.
1 | > summary(pro_model1) |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | Call: lm(formula = lcavol ~ age + lbph + lcp + gleason, data = prostate) Residuals: Min 1Q Median 3Q Max -1.93235 -0.43727 0.00085 0.50103 2.51215 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.78185 1.15927 -0.674 0.502 age 0.02085 0.01315 1.586 0.116 lbph -0.01689 0.06528 -0.259 0.796 lcp 0.51970 0.07400 7.023 3.66e-10 *** gleason 0.13253 0.14733 0.899 0.371 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.8677 on 92 degrees of freedom Multiple R-squared: 0.4805, Adjusted R-squared: 0.458 F-statistic: 21.28 on 4 and 92 DF, p-value: 1.901e-12 |
1 | > summary(pro_model2) |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | Call: lm(formula = lpsa ~ age + lbph + lcp + gleason, data = prostate) Residuals: Min 1Q Median 3Q Max -2.05359 -0.56531 -0.03972 0.57489 2.34593 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.319355 1.277511 1.033 0.3044 age 0.003393 0.014489 0.234 0.8154 lbph 0.133996 0.071935 1.863 0.0657 . lcp 0.412253 0.081544 5.056 2.17e-06 *** gleason 0.148513 0.162362 0.915 0.3627 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9563 on 92 degrees of freedom Multiple R-squared: 0.3423, Adjusted R-squared: 0.3137 F-statistic: 11.97 on 4 and 92 DF, p-value: 7.115e-08 |
How can we improve model performance?
The key difference between regression modeling and other machine learning approaches is that regression typically leaves feature selection and model specification to the user.
Adding non-linear relationships
The below will add age and age2 to the lm() formula using expenses ~ age + age2 form. This will allow the model to separate the linear and non-linear impact of age on lcacol and lpsa.
Converting a numeric variable to a binary indicator
A feature called ifelse() will be used to convert a numeric variable to a binary indicator.
Putting it all together
From the above model improvements, we can try to view an improvement in the regression formula: 1. Adding a non-linear term for age 2. Created gleason index of 6
1 2 3 | > prostate$age2 <- prostate$age^2 > prostate$gle13 <- ifelse(prostate$gleason == 6, 1, 0) |
LCAVOL
1 2 3 | > pro_model1l <- lm(lcavol ~ age + age2 + lbph + lcp + gleason + gle13, data = prostate) > summary(pro_model1l) |