Tractors & Wine Orchards Data Analysis; April 6, 2024
John A. Oludare
2024-04-14
Load the necessary packages
library(readr)
library(ggplot2)Read the dataset
library(readr)
Machinery_STA4813 <- read_csv("~/Documents/Unisa 2024/STA4813/Assignment 1 - 2024/Machinery_STA4813.csv")## Rows: 27 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Type_Tractor, Type_Fruit
## dbl (11): Tractor_Power, Average_Purchase_Price, Salvage_Value, Average_Inve...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Machinery_STA4813$Type_Tractor <- as.factor(Machinery_STA4813$Type_Tractor)
Machinery_STA4813$Type_Fruit <- as.factor(Machinery_STA4813$Type_Fruit)
Machinery_STA4813## # A tibble: 27 × 13
## Tractor_Power Type_Tractor Type_Fruit Average_Purchase_Price Salvage_Value
## <dbl> <fct> <fct> <dbl> <dbl>
## 1 48 2 WD Vineyard 494800 49480
## 2 36 2 WD Orchards 344400 34440
## 3 44 2 WD Orchards 433750 43375
## 4 48 2 WD Orchards 412902 41290
## 5 51 2 WD Orchards 573500 57350
## 6 52 2 WD Orchards 350625 35063
## 7 53 2 WD Orchards 460200 46020
## 8 56 2 WD Orchards 847768 84777
## 9 57 2 WD Orchards 604950 60495
## 10 59 2 WD Orchards 478500 47850
## # ℹ 17 more rows
## # ℹ 8 more variables: Average_Investment <dbl>, Depreciation_Costs <dbl>,
## # Insurance_Licence_Costs <dbl>, Interest_Costs <dbl>,
## # Total_Fixed_Costs <dbl>, Repair_Maintenance_Costs <dbl>, Fuel_Costs <dbl>,
## # Fuel_Usage <dbl>The dataset “Machinery.csv” has undergone cleaning and pre-processing using SPSS and R. It has now been successfully loaded into R-studio using the provided codes. A cursory examination of the data reveals that there are 27 observations and 13 variables, which will be utilized in later analysis.
Using the following R-codes, we will now examine the dataset more closely, including the characteristics of its variables, to comprehend its nature and derive useful insights from our anticipated data analyses.
Explore the dataset
head(Machinery_STA4813)## # A tibble: 6 × 13
## Tractor_Power Type_Tractor Type_Fruit Average_Purchase_Price Salvage_Value
## <dbl> <fct> <fct> <dbl> <dbl>
## 1 48 2 WD Vineyard 494800 49480
## 2 36 2 WD Orchards 344400 34440
## 3 44 2 WD Orchards 433750 43375
## 4 48 2 WD Orchards 412902 41290
## 5 51 2 WD Orchards 573500 57350
## 6 52 2 WD Orchards 350625 35063
## # ℹ 8 more variables: Average_Investment <dbl>, Depreciation_Costs <dbl>,
## # Insurance_Licence_Costs <dbl>, Interest_Costs <dbl>,
## # Total_Fixed_Costs <dbl>, Repair_Maintenance_Costs <dbl>, Fuel_Costs <dbl>,
## # Fuel_Usage <dbl>str(Machinery_STA4813)## spc_tbl_ [27 × 13] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Tractor_Power : num [1:27] 48 36 44 48 51 52 53 56 57 59 ...
## $ Type_Tractor : Factor w/ 2 levels "2 WD","4 WD": 1 1 1 1 1 1 1 1 1 1 ...
## $ Type_Fruit : Factor w/ 2 levels "Orchards","Vineyard": 2 1 1 1 1 1 1 1 1 1 ...
## $ Average_Purchase_Price : num [1:27] 494800 344400 433750 412902 573500 ...
## $ Salvage_Value : num [1:27] 49480 34440 43375 41290 57350 ...
## $ Average_Investment : num [1:27] 222660 154980 195188 185806 258075 ...
## $ Depreciation_Costs : num [1:27] 44.5 31 39 37.2 51.6 ...
## $ Insurance_Licence_Costs : num [1:27] 2.23 1.55 1.95 1.86 2.58 1.58 2.07 3.81 2.72 2.15 ...
## $ Interest_Costs : num [1:27] 28.9 20.1 25.4 24.1 33.5 ...
## $ Total_Fixed_Costs : num [1:27] 46.8 32.5 41 39 54.2 ...
## $ Repair_Maintenance_Costs: num [1:27] 59.4 41.3 52 49.5 68.8 ...
## $ Fuel_Costs : num [1:27] 182 136 166 182 193 ...
## $ Fuel_Usage : num [1:27] 8.64 6.48 7.92 8.64 9.18 ...
## - attr(*, "spec")=
## .. cols(
## .. Tractor_Power = col_double(),
## .. Type_Tractor = col_character(),
## .. Type_Fruit = col_character(),
## .. Average_Purchase_Price = col_double(),
## .. Salvage_Value = col_double(),
## .. Average_Investment = col_double(),
## .. Depreciation_Costs = col_double(),
## .. Insurance_Licence_Costs = col_double(),
## .. Interest_Costs = col_double(),
## .. Total_Fixed_Costs = col_double(),
## .. Repair_Maintenance_Costs = col_double(),
## .. Fuel_Costs = col_double(),
## .. Fuel_Usage = col_double()
## .. )
## - attr(*, "problems")=<externalptr>The “head” command comfortably offers us a preview of the entire data. The variables in this dataset (in order) are: Tractor Power in Kilowatts (kW), Type of Tractor (2-Wheel Drive Tractor or 4-Wheel Drive Tractor), Type of Fruit (Vineyard or Orchards), Average Purchase Price in South African Rand (R), Salvage Value in South African Rand, Average investment in South African Rand, Depreciation costs per hour (in R/h), Insurance and licence costs per hour (in R/h), Interest cost per hour (in R/h), Total fixed costs excluding interest, Repair & Maintenance costs per hour (in R/h), and Fuel costs (in R/h).
The “str” command is another useful way to examine the dataset as an R object. We can clearly see that 2 of the 13 variables are categorical variables and are therefore “Factor Variables”. These 2 are: “Type of Tractor (2-Wheel Drive Tractor or 4-Wheel Drive Tractor)” and “Type of Fruit (Vineyard or Orchards)”.
The remaining 11 variables are quantivative. Of these quantitative variables, “Tractor Power in Kilowatts (kW)”, “Average Purchase Price in South African Rand (R)”, “Salvage Value in South African Rand”, and “Average investment in South African Rand” are integers, while “Depreciation costs per hour (in R/h)”, “Insurance and licence costs per hour (in R/h)”, “Interest cost per hour (in R/h)”, “Total fixed costs excluding interest”, “Repair & Maintenance costs per hour (in R/h)”, and “Fuel costs (in R/h)” are interval or continuous variables.
Initial Data Analysis
summary(Machinery_STA4813)## Tractor_Power Type_Tractor Type_Fruit Average_Purchase_Price
## Min. :36.00 2 WD:11 Orchards:25 Min. : 266650
## 1st Qu.:48.00 4 WD:16 Vineyard: 2 1st Qu.: 446975
## Median :53.00 Median : 526400
## Mean :53.78 Mean : 566553
## 3rd Qu.:59.00 3rd Qu.: 631840
## Max. :71.00 Max. :1171456
## Salvage_Value Average_Investment Depreciation_Costs Insurance_Licence_Costs
## Min. : 26665 Min. :119993 Min. : 24.00 Min. :1.200
## 1st Qu.: 44698 1st Qu.:201139 1st Qu.: 40.23 1st Qu.:2.010
## Median : 52640 Median :236880 Median : 47.38 Median :2.370
## Mean : 56655 Mean :254949 Mean : 50.99 Mean :2.549
## 3rd Qu.: 63184 3rd Qu.:284328 3rd Qu.: 56.87 3rd Qu.:2.840
## Max. :117146 Max. :527155 Max. :105.43 Max. :5.270
## Interest_Costs Total_Fixed_Costs Repair_Maintenance_Costs Fuel_Costs
## Min. :15.60 Min. : 25.20 Min. : 32.00 Min. :136.1
## 1st Qu.:26.14 1st Qu.: 42.24 1st Qu.: 53.63 1st Qu.:181.5
## Median :30.79 Median : 49.74 Median : 63.17 Median :200.4
## Mean :33.14 Mean : 53.54 Mean : 67.99 Mean :203.4
## 3rd Qu.:36.96 3rd Qu.: 59.70 3rd Qu.: 75.82 3rd Qu.:223.1
## Max. :68.53 Max. :110.70 Max. :140.57 Max. :268.5
## Fuel_Usage
## Min. : 6.480
## 1st Qu.: 8.640
## Median : 9.540
## Mean : 9.733
## 3rd Qu.:10.800
## Max. :12.780The output above provides a basic numerical representation of our data. For both categorical variables, we obtain the number of observations and the corresponding class for each factor.
For the numerical variables, we have eleven summary statistics that provide a rough understanding of the distributions.
Finding a linear model to predict tractor average purchase price is our next step.
Question 1(a): Construct appropriate graphs and calculate appropriate summary statistics of both Average_Purchase_Price and Tractor_Power individually and then together. What is your conclusion?
Step 1: Generate graphs for Average_Purchase_Price individually
library(ggplot2)
library(patchwork)
# Create histogram
histogram <- ggplot(Machinery_STA4813, aes(x = Average_Purchase_Price)) +
geom_histogram(binwidth = 100000, fill = "lightblue", color = "black") +
labs(title = "Distribution of Average Purchase Price",
x = "Average Purchase Price",
y = "Frequency")
# Create boxplot
boxplot <- ggplot(Machinery_STA4813, aes(x = "", y = Average_Purchase_Price)) +
geom_boxplot(fill = "lightblue", color = "black") +
labs(title = "Boxplot of Average Purchase Price",
x = "",
y = "Average Purchase Price")
# Create kernel density plot
plot(density(Machinery_STA4813$Average_Purchase_Price),main="Average_Purchase_Price")
rug(Machinery_STA4813$Average_Purchase_Price)
# Arrange the plots side by side
(histogram + boxplot)
summary(Machinery_STA4813$Average_Purchase_Price)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 266650 446975 526400 566553 631840 1171456The observed range of average purchase prices is from R266,650 to R1,171,456, with a mean price of R566,553.
It is clear by looking at the histogram and, more specifically, the kernel density plot, that the distribution is right-skewed, (since the mean value is greater than the median value), and that there are two large outliers in the right tail. It is highly doubtful therefore, that this distribution follows a normal assumption. The appearance of two separate values in the boxplot helps to further clarify the outliers.
# Confirming the presence of outlier(s)
# Identify potential outliers
outliers <- boxplot(Machinery_STA4813$Average_Purchase_Price)$out
# Print the identified outliers
print(outliers)## [1] 911397 1171456The R-codes provided us with outputs that allowed us to more precisely identify the outliers described earlier. These outliers are the “Average Purchase Prices” with values of R911,379 and R1,171,456. Currently, it is uncertain if these outliers are influential points that need to be eliminated or if they demand a transformation of our variable. This will be determined as our analyses go.
Step 2: Generate graphs for Tractor_Power individually
library(ggplot2)
library(patchwork)
# Create histogram
histogram <- ggplot(Machinery_STA4813, aes(x = Tractor_Power)) +
geom_histogram(binwidth = 5, fill = "lightblue", color = "black") +
labs(title = "Distribution of Tractor_Power",
x = "Tractor_Power",
y = "Frequency")
# Create boxplot
boxplot <- ggplot(Machinery_STA4813, aes(x = "", y = Tractor_Power)) +
geom_boxplot(fill = "lightblue", color = "black") +
labs(title = "Boxplot of Tractor_Power",
x = "",
y = "Tractor_Power")
# Create kernel density plot
plot(density(Machinery_STA4813$Tractor_Power),main="Tractor_Power")
rug(Machinery_STA4813$Tractor_Power)
# Arrange the plots side by side
(histogram + boxplot)
summary(Machinery_STA4813$Tractor_Power)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 36.00 48.00 53.00 53.78 59.00 71.00According to the output above, the minimum tractor power is 36.00 kilowatts, and the highest is 71.00 kilowatts. The average tractor power is 53.78 kilowatts. The kernel density figure shows that the data is approximately symmetrically (or uniformly) distributed, with the majority of it centered around the mean. The histogram shows no substantial skewness. The boxplot is also positioned centrally in the data, indicating the absence of outliers.
Step 3: Any meaningful relationships between Average purchase price & tractor power?
# Generate a scatter plot for Average_Purchase_Price and Tractor_Power together
ggplot(Machinery_STA4813, aes(x = Tractor_Power, y = Average_Purchase_Price)) +
geom_point(color = "blue") +
labs(title = "Scatter Plot of Average_Purchase_Price vs. Tractor_Power",
x = "Tractor Power",
y = "Average Purchase Price")
# Calculate correlation coefficient
round(cor(Machinery_STA4813$Tractor_Power, Machinery_STA4813$Average_Purchase_Price), 3)## [1] 0.606A visual representation of the relationship between the two variables is depicted above according to the scatter plot. The data points exhibit an upward trend which signifies a positive correlation.This means that higher-powered tractors tend to be purchased for a greater average price.
The correlation coefficient of 0.606 suggests that Average Purchase Price and Tractor Power are correlated in a moderately positive manner, even though the relationship is not a perfectly positive one.
We should bear in mind however, that correlation does not necessarily indicate causation; i.e. although a positive correlation exists between tractor power and average purchase price, there may be additional determinants of average purchase price that are not accounted for in this analysis.
Question 1(b): Fit a simple linear regression model with Average_Purchase_Price as the response variable and Tractor_Power as the explanatory variable (Model 1). Interpret the estimate slope coefficient β1. Are all assumptions of the linear regression of the met by the data?
# Fitting the linear regression model
model1 <- lm(Average_Purchase_Price ~ Tractor_Power, data = Machinery_STA4813)
# Printing the model summary
model1sum <- summary(model1)
model1sum##
## Call:
## lm(formula = Average_Purchase_Price ~ Tractor_Power, data = Machinery_STA4813)
##
## Residuals:
## Min 1Q Median 3Q Max
## -360398 -104736 7289 49869 429788
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -203953 204543 -0.997 0.328262
## Tractor_Power 14328 3760 3.811 0.000804 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 161100 on 25 degrees of freedom
## Multiple R-squared: 0.3675, Adjusted R-squared: 0.3422
## F-statistic: 14.52 on 1 and 25 DF, p-value: 0.0008037# Printing the coefficients of model1
coef(model1)## (Intercept) Tractor_Power
## -203952.56 14327.59The summary of the linear regression model (model1) is provided above.
Interpretation of the slope coefficient[Tractor_Power (β1)]: The predicted slope coefficient of 14328 represents the expected rise in the average purchase price for each unit of tractor power. In this scenario, as tractor power increases by one unit, the average purchase price is predicted to rise by about R14328, while other factors remain constant. The p-value of 0.000804 indicates that the slope coefficient is statistically significant at 0.05 alpha level.
model1sum$r.squared## [1] 0.3674703The coefficient of determination (R-squared) is 0.3675, indicating that 36.75% of the variability in the model can be explained by the Tractor power variable. Given the relatively low R-squared value, it suggests that the model may not comprehensively incorporate all the factors that contribute to the variability in the response variable. The comparatively low R-squared score suggests that there is room for improvement in the model.
deviance(model1)## [1] 648985891470To further assess the goodness of fit of the regression model1, we have provided the deivance; a measure of the discrepancy between the observed data and the model’s predicted values. The deviance above seems too large for a good model.
“Centering” the mean.
Since the intercept (β0) = -203953; a result which has no contextual meaning when estimating the Average Purchase Price of a tractor, we may want to consider the following improvement analysis:
Machinery_STA4813$Average_Purchase_Price_Centering <- Machinery_STA4813$Average_Purchase_Price - mean( Machinery_STA4813$Average_Purchase_Price)
Machinery_STA4813$Tractor_Power_Centering <- Machinery_STA4813$Tractor_Power - mean(Machinery_STA4813$Tractor_Power)
Model_1_centering <- lm(Average_Purchase_Price_Centering ~ Tractor_Power_Centering, Machinery_STA4813)
summary(Model_1_centering)##
## Call:
## lm(formula = Average_Purchase_Price_Centering ~ Tractor_Power_Centering,
## data = Machinery_STA4813)
##
## Residuals:
## Min 1Q Median 3Q Max
## -360398 -104736 7289 49869 429788
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.991e-11 3.101e+04 0.000 1.000000
## Tractor_Power_Centering 1.433e+04 3.760e+03 3.811 0.000804 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 161100 on 25 degrees of freedom
## Multiple R-squared: 0.3675, Adjusted R-squared: 0.3422
## F-statistic: 14.52 on 1 and 25 DF, p-value: 0.0008037coef(summary(Model_1_centering))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.991025e-11 31007.43 -1.287119e-15 1.000000000
## Tractor_Power_Centering 1.432759e+04 3759.52 3.811015e+00 0.000803672By centering the mean as observed previously, the updated intercept (β0) is now -3.991025e-11 (approximately zero), which provides a more significant interpretation to the Average Purchase Price when the tractor power remains constant (or zero), specifically R0.00 South African Rands will be expended when the Tractor Power is also zero kilowatts. Nevertheless, the remaining parameters (“Multiple R-squared”, “β1”, and standard error) of this model remain unchanged. Therefore, we would pursue additional enhancements for this model.
deviance(Model_1_centering)## [1] 648985891470Deviance still remain unchanged also.
# Diagnostic plots.
# Set the layout to have two rows and two columns
par(mfrow = c(2, 2))
# Scatter plot of residuals vs. fitted values
plot(model1, which = 1)
# Scatter plot of square-root of standardized residuals vs. fitted values
plot(model1, which = 3)
# Histogram of residuals
plot(model1, which = 2)
# Plot of standardized residuals vs leverage
plot(model1, which = 5)
# Reset the layout back to the default
par(mfrow = c(1, 1))# Calculating Cook's distance
cooksd <- cooks.distance(model1)
# Identifying influential points
influential_points <- which(cooksd > 4 / length(cooksd))
# Printing the influential points
print(influential_points)## 26
## 26The first plot displays the residuals plotted against the fitted values. This plot can be utilized to identify the presence of lack of fit and to evaluate the assumption of constant variance. Firstly, there are outliers identified and one of them is even an influential point. The curve has a subtle curved pattern or divergence from the dotted line, as the red line is consistently below the zero line and the residuals tend to be negative. This indicates that the model tends to underestimate the response variable, suggesting that a modification to the model may be necessary, such as transforming one of the variables. The residuals do not appear to have a constant variance.
The next diagnostic tool to explore is the normal quantile-quantile (Q-Q) plot. In this instance, the figure demonstrates that the data distribution closely aligns with the dotted lines, except for the presence of outliers or influential points. The assumption of normality is consequently broken (or violated) as well.
Analyzing the spread and dispersion of data points allows for an interpretation of the scale-location plot that is comparable to that of the residual versus fitted plot. When we look at the spread from left to right in the plot above, we see a funnel shape or a somewhat different spread and pattern, which suggests heteroscedasticity (i.e. non-constant variance). Significant deviations from the general pattern, or outlying points, also signal to important findings or possible problems with the model. In conclusion, there is a breach of the constant variance assumption.
The residual versus leverage plot clearly indicates the presence of an influential observation, suggesting that the model should be modified or a variable should be transformed.
To summarize, despite a modest level of linearity (correlation = 0.606) between the average purchase price and tractor power, the criteria of normality and homoscedasticity are violated. As a result, not all requirements for linear regression are met by the data.
Question 1(c):Create a new variable logAverage_Purchase_Price, the natural log of Average_Purchase_Price. Fit Model 2 where logAverage_Purchase_Price is now the response variable and Tractor_Power is still the explanatory variable. Write down the estimated regression line equation.
# Create the new variable logAverage_Purchase_Price
Machinery_STA4813$logAverage_Purchase_Price <- log(Machinery_STA4813$Average_Purchase_Price)
# Fit Model 2 with logAverage_Purchase_Price as the response variable and Tractor_Power as the explanatory variable
model2 <- lm(logAverage_Purchase_Price ~ Tractor_Power, data = Machinery_STA4813)
# Print the estimated regression line equation
summary(model2)##
## Call:
## lm(formula = logAverage_Purchase_Price ~ Tractor_Power, data = Machinery_STA4813)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.80020 -0.13471 0.01824 0.12945 0.49020
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.918847 0.339158 35.142 < 2e-16 ***
## Tractor_Power 0.023708 0.006234 3.803 0.00082 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2672 on 25 degrees of freedom
## Multiple R-squared: 0.3665, Adjusted R-squared: 0.3412
## F-statistic: 14.46 on 1 and 25 DF, p-value: 0.00082coef(model2)## (Intercept) Tractor_Power
## 11.91884732 0.02370773logAverage_Purchase_Price = 11.918847 + 0.023708 * Tractor_Power
In the equation above, the intercept (estimated coefficient for the (Intercept) term) is 11.918847, and the slope (estimated coefficient for the Tractor_Power term) is 0.023708.
deviance(model2)## [1] 1.784321Question 1(d): How does logAverage_Purchase_Price change when Tractor_Power increases by 1?
The coefficient indicates that for every one-unit rise in the Tractor_Power, the logAverage_Purchase_Price is predicted to increase by about 0.023708.
Question 1(e): Recall that log(a) −log(b) = log( a/b). Use this to explain/find out how Average_Purchase_Price changes as Tractor_Power increases by 1.
Since logAverage_Purchase_Price = 11.918847 + 0.023708 * Tractor_Power, then as tractor power increase by 1, it implies: log(Average_Purchase_Price_new) = 11.918847 + 0.023708 * 1 = 11.94255. Hence,
Average_Purchase_Price_new <- exp(11.94255)
Average_Purchase_Price_new## [1] 153668Therefore average power changes by 153668 units.
Question 1(f): Are linear least squares regression (LLSR) assumptions satisfied in Model 2? Why or why not?
Just as we diagnosed for model1
# Diagnostic plots.
# Set the layout to have two rows and two columns
par(mfrow = c(2, 2))
# Scatter plot of residuals vs. fitted values
plot(model2, which = 1)
# Scatter plot of square-root of standardized residuals vs. fitted values
plot(model2, which = 3)
# Histogram of residuals
plot(model2, which = 2)
# Plot of standardized residuals vs leverage
plot(model2, which = 5)
# Reset the layout back to the default
par(mfrow = c(1, 1))Again,
The first plot displays the residuals plotted against the fitted values. The curve has a no divergence from the dotted line this time, as the red line is consistently along line and the residuals randomly distributed. This indicates that the model tends to represented the response variable. There might be a worrying outlier identified as 21. The residuals appear to have a constant variance.
The normal quantile-quantile (Q-Q) plot in this case, demonstrates that the data distribution closely aligns with the dotted lines. However, some slight cyclic patterns are observed. The data appears normally distributed.
The scale-location plot looks at the spread from left to right, and we see a more even spread, which suggests homoscedasticity (constant variance). However, there might be a prominent outlier which in conclusion, shows a breach of the constant variance assumption.
The residual versus leverage plot indicates the absence of an influential observation, suggesting that the model satisfies normality.
To summarize, the criterion of normality is satisfied, but that of homoscedasticity seems violated. As a result, not all requirements for linear regression are met by the data.
_** Question 1(g):Fit Model 3 with Average_Purchase_Price as the response and Tractor_Power and Type_Tractor as predictors or explanatory variables. Write down the estimated regression line equation and interpret all the regression coefficients in the fitted model**_
# Convert 2WD to 0 and 4WD to 1
Machinery_STA4813$Type_Tractor <- as.numeric(Machinery_STA4813$Type_Tractor) - 1
Machinery_STA4813$Type_Tractor## [1] 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1model_3 <- lm(Average_Purchase_Price ~ Tractor_Power + Type_Tractor, data = Machinery_STA4813)
summary(model_3)##
## Call:
## lm(formula = Average_Purchase_Price ~ Tractor_Power + Type_Tractor,
## data = Machinery_STA4813)
##
## Residuals:
## Min 1Q Median 3Q Max
## -373071 -110850 22677 38741 421535
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -196061 207907 -0.943 0.35507
## Tractor_Power 13775 3939 3.497 0.00186 **
## Type_Tractor 36836 66127 0.557 0.58265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 163400 on 24 degrees of freedom
## Multiple R-squared: 0.3755, Adjusted R-squared: 0.3235
## F-statistic: 7.217 on 2 and 24 DF, p-value: 0.003516coef(model_3)## (Intercept) Tractor_Power Type_Tractor
## -196060.95 13774.94 36835.91# Extract the coefficients and their variable names
coefficients <- coef(model_3)
variable_names <- names(coefficients)
# Create the MLR equation string
mlr_equation <- paste("Average_Purchase_Price =", coefficients[1], "+", paste(coefficients[-1], variable_names[-1], collapse = " + "))
# Print the MLR equation
cat("MLR Equation:", mlr_equation)## MLR Equation: Average_Purchase_Price = -196060.954444562 + 13774.9388278195 Tractor_Power + 36835.9120005814 Type_Tractormodel_3sum<-summary(model_3)
model_3sum$r.squared## [1] 0.3755442The intercept term (β0 = -159225.04) indicates the projected average purchase price for a tractor that has no Tractor_Power and is of Type_Tractor 2WD. Nevertheless, due to the fact that 2WD is not a numerical variable, we approach the interpretation of this coefficient with care and place greater emphasis on the other predictors.
The coefficient Tractor_Power (β1 = 13774.94) represents the expected rise in the average purchase price for each incremental unit of Tractor_Power, assuming that the Type_Tractor remains unchanged. Hence, a one-unit increment in Tractor_Power is projected to result in an average purchase price increase of around R13,774.94.
Type_Tractor (β2 = 36835.91): This coefficient quantifies the expected variation in the mean buying cost between Type_Tractor 4WD and Type_Tractor 2WD, while keeping Tractor_Power constant. On average, Type_Tractor 4WD tractors are anticipated to have a higher mean purchase price than Type_Tractor 2WD tractors by roughly R36,835.91.
_** Question 1(h): Are LLSR assumptions satisfied in Model 3? Why or why not?**_
# Diagnostic plots.
# Set the layout to have two rows and two columns
par(mfrow = c(2, 2))
# Scatter plot of residuals vs. fitted values
plot(model_3, which = 1)
# Scatter plot of square-root of standardized residuals vs. fitted values
plot(model_3, which = 3)
# Histogram of residuals
plot(model_3, which = 2)
# Plot of standardized residuals vs leverage
plot(model_3, which = 5)
# Reset the layout back to the default
par(mfrow = c(1, 1))_** Question 1(i): Calculate the correlation matrix of Average_Purchase_Price, Fuel_Usage, Salvage_Value, Average_Investment, Depreciation_Costs, Insurance_Licence_Costs, Interest_Costs and Repair_Maintenance_Costs. What do you observe in the correlation matrix? Do you recommend any explanatory variable from these variables (from Fuel_Usage, Salvage_Value, Average_Investment, Depreciation_Costs, Insurance_Licence_Costs, Interest_Costs and Repair_Maintenance_Costs) to be added to Model 3 to create a model with a higher adjusted R2 value? Why? Explain the reason from LLSR assumptions point of view**_
correlation_matrix <- cor(Machinery_STA4813[c("Average_Purchase_Price", "Fuel_Usage", "Salvage_Value", "Average_Investment", "Depreciation_Costs", "Insurance_Licence_Costs", "Interest_Costs", "Repair_Maintenance_Costs")])
correlation_matrix## Average_Purchase_Price Fuel_Usage Salvage_Value
## Average_Purchase_Price 1.0000000 0.6524046 1.0000000
## Fuel_Usage 0.6524046 1.0000000 0.6524043
## Salvage_Value 1.0000000 0.6524043 1.0000000
## Average_Investment 1.0000000 0.6524046 1.0000000
## Depreciation_Costs 1.0000000 0.6524095 1.0000000
## Insurance_Licence_Costs 0.9999970 0.6517110 0.9999970
## Interest_Costs 1.0000000 0.6523766 1.0000000
## Repair_Maintenance_Costs 1.0000000 0.6524099 1.0000000
## Average_Investment Depreciation_Costs
## Average_Purchase_Price 1.0000000 1.0000000
## Fuel_Usage 0.6524046 0.6524095
## Salvage_Value 1.0000000 1.0000000
## Average_Investment 1.0000000 1.0000000
## Depreciation_Costs 1.0000000 1.0000000
## Insurance_Licence_Costs 0.9999970 0.9999969
## Interest_Costs 1.0000000 1.0000000
## Repair_Maintenance_Costs 1.0000000 1.0000000
## Insurance_Licence_Costs Interest_Costs
## Average_Purchase_Price 0.9999970 1.0000000
## Fuel_Usage 0.6517110 0.6523766
## Salvage_Value 0.9999970 1.0000000
## Average_Investment 0.9999970 1.0000000
## Depreciation_Costs 0.9999969 1.0000000
## Insurance_Licence_Costs 1.0000000 0.9999971
## Interest_Costs 0.9999971 1.0000000
## Repair_Maintenance_Costs 0.9999971 1.0000000
## Repair_Maintenance_Costs
## Average_Purchase_Price 1.0000000
## Fuel_Usage 0.6524099
## Salvage_Value 1.0000000
## Average_Investment 1.0000000
## Depreciation_Costs 1.0000000
## Insurance_Licence_Costs 0.9999971
## Interest_Costs 1.0000000
## Repair_Maintenance_Costs 1.0000000# Calculate the absolute correlation values
abs_correlation <- abs(correlation_matrix[, "Average_Purchase_Price"])
# Find the variable with the least correlation
variable_least_correlation <- names(abs_correlation)[which.min(abs_correlation)]
# Print the variable with the least correlation
cat("Variable with the least correlation:", variable_least_correlation)## Variable with the least correlation: Fuel_UsageFuel_Usage has the least correlation and can be added to model_3. Almost every variable risk multicollinearity because of their very high or perfect correlation with one another. Hence,
model_3_new <- lm(Average_Purchase_Price ~ Tractor_Power + Type_Tractor + Fuel_Usage, data = Machinery_STA4813)
summary(model_3_new)##
## Call:
## lm(formula = Average_Purchase_Price ~ Tractor_Power + Type_Tractor +
## Fuel_Usage, data = Machinery_STA4813)
##
## Residuals:
## Min 1Q Median 3Q Max
## -351183 -88435 17233 54730 451015
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -175796 180511 -0.974 0.34025
## Tractor_Power -65239 26730 -2.441 0.02278 *
## Type_Tractor 49718 57535 0.864 0.39643
## Fuel_Usage 433696 145515 2.980 0.00669 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 141800 on 23 degrees of freedom
## Multiple R-squared: 0.5495, Adjusted R-squared: 0.4908
## F-statistic: 9.352 on 3 and 23 DF, p-value: 0.0003143model_3_newsum<-summary(model_3_new)
model_3_newsum$adj.r.squared## [1] 0.4907664There is a rise in the new adjusted R-squared value after adding fuel usage.
# Diagnostic plots.
# Set the layout to have two rows and two columns
par(mfrow = c(2, 2))
# Scatter plot of residuals vs. fitted values
plot(model_3_new, which = 1)
# Scatter plot of square-root of standardized residuals vs. fitted values
plot(model_3_new, which = 3)
# Histogram of residuals
plot(model_3_new, which = 2)
# Plot of standardized residuals vs leverage
plot(model_3_new, which = 5)
# Reset the layout back to the default
par(mfrow = c(1, 1))The LLSR assumptions are still not met primarily because of the presence of the outliers or influential point.
We will now check for the outlier and remove it.
# Assuming you have previously fitted the regression model "model_3_new"
# Calculating Cook's distance
cooksd <- cooks.distance(model_3_new)
# Identifying influential points
influential_points <- which(cooksd > 4 / length(cooksd))
# Identifying outliers (including influential points)
outliers <- which(cooksd > 4 / length(cooksd))
# Printing the influential points
print(influential_points)## 26
## 26# Printing the outliers
print(outliers)## 26
## 26Removal of influential observation
# Assuming the outlier(s) or influential observation(s) are identified and stored in the "influential_points" variable
# Create a new dataset without the outlier(s)
Machinery_STA4813_New <- Machinery_STA4813[!rownames(Machinery_STA4813) %in% influential_points, ]
# Re-fit the regression model with the new dataset
model4 <- lm(Average_Purchase_Price ~ Tractor_Power + Type_Tractor + Fuel_Usage, data = Machinery_STA4813_New)
# Print the summary of the new model
summary(model4)##
## Call:
## lm(formula = Average_Purchase_Price ~ Tractor_Power + Type_Tractor +
## Fuel_Usage, data = Machinery_STA4813_New)
##
## Residuals:
## Min 1Q Median 3Q Max
## -308927 -35590 2832 54140 143448
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11540 134156 -0.086 0.932232
## Tractor_Power -75185 19310 -3.894 0.000781 ***
## Type_Tractor 31767 41490 0.766 0.452012
## Fuel_Usage 470888 104790 4.494 0.000181 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 101800 on 22 degrees of freedom
## Multiple R-squared: 0.6471, Adjusted R-squared: 0.599
## F-statistic: 13.45 on 3 and 22 DF, p-value: 3.364e-05model4sum<-summary(model4)
model4sum$adj.r.squared## [1] 0.5989653The model now looks better with an improved Adjusted R-squared value of approximately 60%.
# Diagnostic plots.
# Set the layout to have two rows and two columns
par(mfrow = c(2, 2))
# Scatter plot of residuals vs. fitted values
plot(model4, which = 1)
# Scatter plot of square-root of standardized residuals vs. fitted values
plot(model4, which = 3)
# Histogram of residuals
plot(model4, which = 2)
# Plot of standardized residuals vs leverage
plot(model4, which = 5)
# Reset the layout back to the default
par(mfrow = c(1, 1))The LLSR assumptions are now fulfilled as seen from the four diagnostic plots above; Linearity, Normality, Homoscedasticity and independence of errors are all valid.