Labour Estimating & Performance Review
2024–2025 Patio Installations | Statistical Methodology
October 27, 2026
Methodology at a Glance
This document outlines the statistical methodology supporting the companion Labour Estimating & Performance Review management decision brief.
The analysis examines 300 patio installations across four crews during the 2024 and 2025 construction seasons, using simulated operational records supplemented with historical Ottawa weather data.
Our methodology consisted of four stages:
- Exploratory Data Analysis: Examined company-wide and crew-level labour overruns to identify patterns and differences in estimating accuracy.
- Regression Analysis: Used multiple linear regression to investigate associations between labour estimating discrepancies, crew assignment, project characteristics, and operational complications.
- Model Diagnostics: Evaluated regression assumptions, statistical significance, and influential observations through diagnostic tests and HC3 robust inference.
- Sensitivity Analysis: Assessed the stability of the principal findings by examining how selected influential projects affected the estimated regression coefficients.
Analytical objective: Identify factors associated with labour estimating discrepancies, assess the reliability of those associations, and highlight projects warranting management investigation.
The model was developed primarily for explanation and decision support rather than individual project-level prediction.
Data Acquisition and Synthetic Data Generation
Overview of the Dataset
The primary analytical dataset, patio_projects_synthetic.csv, was generated in R using a purpose-built simulation program.
The dataset contains 300 fictional interlock patio installations, divided equally between the 2024 and 2025 landscaping seasons, with 150 projects assigned to each year.
Projects were allocated among four crews, labelled A, B, C, and D. Each crew was assumed to consist of four workers.
The unit of observation is an individual patio installation. Each record contains information about project dimensions, site accessibility, construction requirements, crew assignment, scheduling, weather exposure, recorded operational complications, and estimated and actual labour requirements.
The data-generation process was designed to resemble several features of residential interlock installation work without reproducing any real company’s operational records.
The simulation used a fixed random seed (20261021) to make the randomly generated observations reproducible when the same input files, estimation functions, and software environment are used.
Generation of Project Characteristics
Project characteristics were generated before calculating actual labour requirements.
Patio areas were randomly generated between approximately 200 and 1,000 square feet, using a log-uniform distribution. This provided installations of different sizes rather than assuming all projects were concentrated around a single patio area.
Each patio was assigned a simulated length and width based on its area and a randomly generated aspect ratio.
Additional construction characteristics were generated to reflect common differences between residential interlock projects.
These included whether the installation was constructed against a house, whether existing interlock had to be removed, whether the site permitted skid-steer access, the method used for removing existing pavers, the type of lawn restoration required, and the complexity of cutting and fitting pavers.
Several generation probabilities were specified in advance. For example, approximately 28% of projects were assigned replacement work, 85% were assumed to be constructed against a house, and 32% were assigned sod rather than grass-seed restoration.
Skid-steer accessibility was generated using a probability that varied with patio area and replacement status. This allowed some project characteristics to be related rather than treating every variable as independent.
Cutting complexity was also allowed to vary across installations, using a positive random multiplier restricted to a range of 0.90 to 1.65.
Crew Assignment
Projects were assigned to one of four fictional installation crews.
Assignment was randomized but not completely independent of project characteristics. This was intended to represent a company in which crews perform overlapping types of work but may be somewhat more likely to receive particular assignments.
For example, the simulation made Crew B somewhat more likely to receive projects with restricted skid-steer access, Crew C somewhat more likely to receive patios exceeding 650 square feet, and Crew D somewhat more likely to receive replacement installations.
These differences were deliberately modest, allowing each crew to complete a variety of projects.
This feature was important for the subsequent statistical analysis. Directly comparing crews without accounting for differences in their assigned work could produce misleading conclusions about their labour requirements.
Construction of Original Labour Estimates
Rather than randomly assigning estimated labour hours, we developed a separate, rule-based labour-estimating function in R, estimate_patio(), contained in labor_estimation.R.
The function was constructed around a reference installation consisting of a 400-square-foot patio (20 × 20 feet), built against a house by a four-person crew, with skid-steer access and grass-seed restoration.
For this reference installation, the labour-estimating model allocated approximately 11.71 active crew-elapsed hours, equivalent to approximately 46.83 person-hours.
The installation process was divided into ten construction stages:
- Mobilization and layout
- Existing paver removal, where applicable
- Excavation and hauling
- Granular A placement and compaction
- Bedding preparation and grading
- Transporting and laying full pavers
- Cutting and fitting
- Edge-restraint installation
- Polymeric sand application and final compaction
- Lawn restoration
The estimated duration of each stage was adjusted according to relevant project characteristics.
For example, excavation requirements varied with the area requiring preparation, while cutting requirements depended partly on patio perimeter and cutting complexity. Restricted equipment access increased the assumed labour requirements for certain material-handling and installation activities. Replacement projects also required additional removal work.
Some activities included a fixed time component, while others varied primarily with project size or geometric characteristics.
The project-level estimate was calculated by summing the estimated durations of all applicable construction stages.
Labour hours were defined as active on-site crew-elapsed hours, excluding scheduled breaks, travel from the company’s yard, and the separate truck driver’s time.
Because all simulated crews contained four workers, total person-hours were calculated by multiplying active crew-elapsed hours by four.
The original estimates were calculated separately from actual labour outcomes. This distinction allowed the simulation to represent an estimating process in which contractors make decisions using information available before construction, while actual labour requirements may subsequently differ.
Project Scheduling and Historical Weather Data
The simulation incorporated historical daily weather observations from Environment and Climate Change Canada’s Historical Climate Data service, using Ottawa International Airport (station ID 49568) for 2024 and 2025 (Environment and Climate Change Canada, n.d.).
The downloaded observations were processed in R and stored in ottawa_weather_clean.csv. The dataset included daily precipitation and maximum temperature, which were used to characterize weather exposure during the simulated construction schedules.
The weather dataset included daily precipitation and maximum temperature. These variables were used to characterize weather exposure during the fictional construction schedule.
Projects were assigned planned workdays between April and November, excluding weekends.
For scheduling purposes, the simulation assumed approximately 9.1 active crew hours per workday. The number of planned workdays was determined from each project’s original labour estimate.
The scheduling algorithm also prevented a crew from being assigned to two planned projects on the same workday.
Each project’s planned dates were then matched with the corresponding Ottawa weather observations.
Two weather indicators were constructed:
- Wet-weather fraction: The proportion of planned workdays with at least 5 mm of precipitation.
- Hot-weather fraction: The proportion of planned workdays with a maximum temperature of at least 30°C.
For example, a project with two qualifying rainy days across four planned workdays received a wet-weather fraction of 0.50.
A small number of missing precipitation or maximum-temperature observations were replaced using the median observed value for the corresponding year and month. Dates requiring imputation were flagged in the data.
Weather exposure was based on planned workdays rather than reconstructed actual construction dates. The simulation therefore does not establish the precise weather conditions that a real crew would have experienced after unexpected schedule changes.
Generation of Operational Complications
Operational complications were introduced to represent circumstances that could increase labour requirements after an estimate had been prepared.
Three principal indicators were generated:
- Difficult soil: Assigned with a probability of 18%.
- Utility conflicts: Assigned with a probability of 12% for replacement projects and 7% for new installations.
- Delivery issues: Assigned with a probability of 7%.
These variables were generated independently of the final labour outcomes, although some probabilities depended on project characteristics.
The complications were then incorporated into the process used to generate actual labour hours.
Their effects were allowed to differ by construction stage. For example, difficult soil and utility conflicts primarily affected excavation, while delivery issues affected selected material-handling and installation stages.
This stage-specific approach was intended to reflect the fact that a disruption may affect some construction activities more than others.
The simulation also allowed occasional unexpected delays beyond the recorded complications. These additional delays were deliberately excluded from the main analytical predictors, introducing situations in which recorded project characteristics would not fully explain actual labour requirements.
Simulation of Actual Labour Requirements
Actual labour hours were generated separately from the original estimates.
For each construction stage, the estimated duration was multiplied by a factor incorporating crew-specific differences, selected operational conditions, weather exposure, and random variation.
The simulation assigned modest underlying crew differences in its data-generating process. Relative to Crew A, the crew-specific log multipliers were approximately -0.05 for Crew B, +0.02 for Crew C, and +0.06 for Crew D. Their contributions varied somewhat across construction stages.
These values were specified during data generation; they are not regression estimates or observed measurements of actual crew productivity.
Operational conditions were also introduced through stage-specific multipliers. For example, difficult soil increased simulated excavation requirements, while wetter conditions affected several stages, including excavation, base preparation, and finishing.
Random variation was introduced at two levels.
First, each project received a common random disturbance that influenced multiple construction stages. Second, additional variation was introduced separately at each stage.
This allowed construction-stage outcomes within the same project to be related while still exhibiting individual variation.
Finally, approximately 5% of projects were assigned an additional unplanned-delay event. When selected, the event added between 1.5 and 4.5 active crew-elapsed hours to one of three construction stages: excavation, granular-base preparation, or paver installation.
Actual project labour was obtained by summing the simulated durations across all construction stages.
Importantly, actual labour hours were not generated directly from the multiple linear regression model used later in the analysis.
Instead, the regression model was fitted to outcomes generated through this more detailed construction-stage process. This allowed the fitted model to be imperfect and left room for residual variation, unmeasured disturbances, and model-specification limitations.
Construction of the Analytical Variables
After generating the estimated and actual labour hours, several analytical variables were calculated.
For project i, let:
- Aᵢ represent actual active crew-elapsed hours.
- Eᵢ represent originally estimated active crew-elapsed hours.
The actual-to-estimated labour ratio was defined as Aᵢ / Eᵢ.
A value greater than 1 indicates that the project exceeded its original labour estimate. A value below 1 indicates that the project required fewer hours than estimated.
The percentage labour overrun was calculated as:
100 × (Aᵢ / Eᵢ − 1)
For the regression analysis, the response variable was defined as the natural logarithm of the actual-to-estimated labour ratio:
log_hours_ratio = log(Aᵢ / Eᵢ)
Using this response allowed the regression to examine proportional differences in labour requirements relative to the original estimates, rather than absolute differences in hours.
Generated Files and Data Validation
The simulation produced three principal output files:
patio_projects_synthetic.csv: Project-level observations used for descriptive analysis and regression modelling.patio_project_stages_synthetic.csv: Estimated and simulated actual labour requirements for individual construction stages.SIMULATION_TRUTH_DO_NOT_USE_AS_PREDICTORS.csv: Internal simulation information, including predefined crew effects and randomly introduced disturbances.
The third file was kept separate from the analytical dataset because it contains information that would not ordinarily be available to management in a real operational investigation.
Using these hidden generating parameters as explanatory variables would provide the analyst with information unavailable in a genuine business setting.
Several validation checks were included in the generation script. These verified that the dataset contained 300 unique projects, that estimated and actual labour values were positive and finite, that the appropriate number of construction-stage records had been generated, and that planned crew assignments did not overlap.
The fixed random seed supported reproducibility of the generated records, subject to using the same input data and simulation-code versions.
Limitations of the Synthetic Data
The synthetic dataset provides a controlled setting for demonstrating an applied statistical investigation, but its limitations must be acknowledged.
The distributions governing project characteristics, crew assignments, operational complications, and labour outcomes were selected as modelling assumptions rather than estimated from historical records supplied by a landscaping company.
Consequently, the simulated frequency and magnitude of labour overruns should not be interpreted as representative of actual landscaping contractors.
In addition, several potentially important operational variables were not explicitly measured, including foreman experience, individual worker experience, changes in crew composition, and detailed task-level interruptions.
Although the simulation incorporates real weather observations, the project schedules and labour outcomes remain fictional.
The exercise is therefore best understood as a methodological demonstration of how a contractor’s project records could be analyzed, rather than evidence that the same patterns or estimated relationships would occur in an actual business.
Exploratory Data Analysis and Model Development
Before fitting the regression models, we conducted an exploratory analysis of estimated and actual labour requirements across the 2024 and 2025 seasons. The descriptive results, including company-wide and crew-level labour overruns, are presented in the companion management decision brief.
The exploratory analysis revealed differences in labour overruns between crews. However, examining project characteristics also indicated that crews were not necessarily assigned equivalent types of work. For example, differences were observed in patio sizes, equipment accessibility, and the frequency of recorded operational complications.
Consequently, direct comparisons of aggregate labour overruns could be misleading, as the observed differences might partly reflect variations in project characteristics rather than differences attributable to crew assignment.
These observations motivated the use of multiple linear regression model.
Regression Specification and Estimation
Response variable
For project \(i\), let \(E_i>0\) denote the original estimated active crew-elapsed hours and \(A_i>0\) the simulated actual active crew-elapsed hours. The response variable is
\[ Y_i=\log\left(\frac{A_i}{E_i}\right). \]
A value of \(Y_i=0\) indicates that actual hours equalled the original estimate; positive values correspond to overruns, and negative values to work completed below the estimate. This response focuses on relative estimating discrepancies instead of absolute hours, which tend to be larger for larger projects simply because those projects require more work.
The logarithm also permits a multiplicative interpretation of regression coefficients. For a binary indicator or a comparison between crews, a coefficient \(\beta\) corresponds to a fitted percentage difference in the actual-to-estimated ratio of
\[ 100\left(e^\beta-1\right)\%. \]
For a continuous predictor changing by \(\Delta\), the corresponding difference is \(100(e^{\beta\Delta}-1)\%\).
Model
Regression Model Specification
Following the exploratory analysis, we fitted a multiple linear regression model to investigate whether differences between crews persisted after accounting for measured project characteristics and operational complications.
The response variable was defined as the natural logarithm of the ratio of actual to originally estimated active crew-elapsed hours:
\[ Y_i=\log\left(\frac{A_i}{E_i}\right), \]
where \(A_i\) and \(E_i\) represent the actual and estimated active crew-elapsed hours, respectively, for project \(i\).
A positive value of \(Y_i\) indicates that the project exceeded its original estimate, while a negative value indicates that fewer hours were required than estimated.
Initial Regression Model: Crew A as the Reference
The initial model included crew assignment, year, project characteristics, recorded operational complications, and weather exposure as explanatory variables.
Crew A was initially used as the reference category.
The regression model was specified as follows:
\[ \begin{aligned} Y_i ={}& \beta_0 +\beta_1 I(\text{Crew}_i=B) +\beta_2 I(\text{Crew}_i=C)\\ &+\beta_3 I(\text{Crew}_i=D) +\beta_4 I(\text{Year}_i=2025)\\ &+\beta_5\log(\text{Area}_i) +\beta_6\text{SkidAccess}_i\\ &+\beta_7\text{Replacement}_i +\beta_8\text{CutComplexity}_i\\ &+\beta_9\text{AgainstHouse}_i +\beta_{10}\text{Sod}_i\\ &+\beta_{11}\text{DifficultSoil}_i +\beta_{12}\text{UtilityConflict}_i\\ &+\beta_{13}\text{DeliveryIssue}_i +\beta_{14}\text{WetFraction}_i\\ &+\beta_{15}\text{HotFraction}_i +\varepsilon_i. \end{aligned} \]
Here, \(I(\cdot)\) denotes an indicator variable equal to one when the specified condition is satisfied and zero otherwise. The remaining binary variables are coded in the same manner. The error term \(\varepsilon_i\) represents the difference between the observed response and its conditional mean under the specified model.
Since Crew A is the reference category, no separate indicator is included for Crew A. The coefficients associated with Crews B, C, and D therefore represent differences in the expected log actual-to-estimated labour ratio relative to Crew A, holding the other explanatory variables constant.
The model was estimated using ordinary least squares (OLS) on all 300 simulated patio installations.
Initial OLS results
| Variable | Estimate | Std. Error | t-value | p-value |
|---|---|---|---|---|
| Intercept | -0.1197692 | 0.0675332 | -1.773 | 0.077220 |
| Crew B | -0.0476194 | 0.0125861 | -3.784 | 0.000189 |
| Crew C | 0.0144514 | 0.0136040 | 1.062 | 0.289006 |
| Crew D | 0.0533372 | 0.0128246 | 4.159 | 0.0000424 |
| Year 2025 | -0.0026335 | 0.0089806 | -0.293 | 0.769549 |
| Log(Patio area) | 0.0156064 | 0.0098782 | 1.580 | 0.115246 |
| Skid-steer access | 0.0062258 | 0.0098947 | 0.629 | 0.529722 |
| Replacement | -0.0007949 | 0.0106141 | -0.075 | 0.940353 |
| Cutting complexity | 0.0293420 | 0.0266263 | 1.102 | 0.271398 |
| Against house | -0.0047350 | 0.0123535 | -0.383 | 0.701792 |
| Sod restoration | -0.0089778 | 0.0094722 | -0.948 | 0.344035 |
| Difficult soil | 0.0432639 | 0.0117135 | 3.693 | 0.000265 |
| Utility conflict | 0.0392998 | 0.0168812 | 2.328 | 0.020614 |
| Delivery issue | 0.0520641 | 0.0169820 | 3.066 | 0.002380 |
| Wet-weather fraction | 0.0932763 | 0.0190266 | 4.902 | 0.00000159 |
| Hot-weather fraction | 0.0132196 | 0.0188692 | 0.701 | 0.484130 |
Model summary
| Statistic | Value |
|---|---|
| Observations | 300 |
| Residual degrees of freedom | 284 |
| Residual standard error | 0.07679 |
| Multiple \(R^2\) | 0.3075 |
| Adjusted \(R^2\) | 0.2709 |
| Overall F-statistic | 8.406 |
| F-statistic degrees of freedom | (15, 284) |
| Overall F-test p-value | \(7.08\times10^{-16}\) |
Note: The standard errors, t-statistics, and p-values in this table are from the conventional OLS summary. Heteroskedasticity-robust inference is examined separately in the regression diagnostics section.
The estimated model explained approximately 30.8% of the observed variation in the logarithm of the actual-to-estimated labour ratio.
Using Crew A as the reference, Crew B had a negative estimated coefficient, while Crews C and D had positive coefficients. This parameterization therefore provided direct comparisons between Crew A and each of the other three crews.
Reparameterization: Crew B as the Reference
Following the initial regression analysis, we changed the reference category from Crew A to Crew B.
This decision was motivated by the exploratory findings. Crew B’s aggregate actual labour requirements remained relatively close to its original estimates in both construction seasons, recording an overrun of approximately 0.20% in 2024 and finishing approximately 1.00% below its estimates in 2025.
Consequently, Crew B was considered a useful operational benchmark for interpreting differences in estimating accuracy.
The purpose of selecting Crew B was not to establish that it was the most productive crew. Rather, it allowed us to compare the other crews against the crew whose aggregate labour requirements were closest to its original estimates.
The regression was therefore reparameterized as:
\[ \begin{aligned} Y_i ={}& \alpha_0 +\alpha_1 I(\text{Crew}_i=A) +\alpha_2 I(\text{Crew}_i=C)\\ &+\alpha_3 I(\text{Crew}_i=D) +\alpha_4 I(\text{Year}_i=2025)\\ &+\alpha_5\log(\text{Area}_i) +\alpha_6\text{SkidAccess}_i\\ &+\alpha_7\text{Replacement}_i +\alpha_8\text{CutComplexity}_i\\ &+\alpha_9\text{AgainstHouse}_i +\alpha_{10}\text{Sod}_i\\ &+\alpha_{11}\text{DifficultSoil}_i +\alpha_{12}\text{UtilityConflict}_i\\ &+\alpha_{13}\text{DeliveryIssue}_i +\alpha_{14}\text{WetFraction}_i\\ &+\alpha_{15}\text{HotFraction}_i +\varepsilon_i. \end{aligned} \]
The explanatory variables are identical to those in the initial model. The only difference is the reference category used for the crew indicators.
Since Crew B is now the reference, the coefficients associated with Crews A, C, and D represent differences relative to Crew B.
OLS results with Crew B as the reference
| Variable | Estimate | Std. Error | t-value | p-value |
|---|---|---|---|---|
| Intercept | -0.1673885 | 0.0676698 | -2.474 | 0.013961 |
| Crew A | 0.0476194 | 0.0125861 | 3.784 | 0.000189 |
| Crew C | 0.0620708 | 0.0128339 | 4.836 | 0.00000217 |
| Crew D | 0.1009566 | 0.0119617 | 8.440 | \(1.65\times10^{-15}\) |
| Year 2025 | -0.0026335 | 0.0089806 | -0.293 | 0.769549 |
| Log(Patio area) | 0.0156064 | 0.0098782 | 1.580 | 0.115246 |
| Skid-steer access | 0.0062258 | 0.0098947 | 0.629 | 0.529722 |
| Replacement | -0.0007949 | 0.0106141 | -0.075 | 0.940353 |
| Cutting complexity | 0.0293420 | 0.0266263 | 1.102 | 0.271398 |
| Against house | -0.0047350 | 0.0123535 | -0.383 | 0.701792 |
| Sod restoration | -0.0089778 | 0.0094722 | -0.948 | 0.344035 |
| Difficult soil | 0.0432639 | 0.0117135 | 3.693 | 0.000265 |
| Utility conflict | 0.0392998 | 0.0168812 | 2.328 | 0.020614 |
| Delivery issue | 0.0520641 | 0.0169820 | 3.066 | 0.002380 |
| Wet-weather fraction | 0.0932763 | 0.0190266 | 4.902 | 0.00000159 |
| Hot-weather fraction | 0.0132196 | 0.0188692 | 0.701 | 0.484130 |
Model summary
| Statistic | Value |
|---|---|
| Observations | 300 |
| Residual degrees of freedom | 284 |
| Residual standard error | 0.07679 |
| Multiple \(R^2\) | 0.3075 |
| Adjusted \(R^2\) | 0.2709 |
| Overall F-statistic | 8.406 |
| F-statistic degrees of freedom | (15, 284) |
| Overall F-test p-value | \(7.08\times10^{-16}\) |
Note: These are conventional OLS coefficient tests. Robust standard errors and their implications for inference are discussed in the subsequent diagnostic analysis.
Equivalence of the Two Parameterizations
Although the coefficient tables differ, the two specifications represent the same fitted regression model.
Changing the reference category changes the interpretation of the intercept and the crew coefficients, but it does not change the estimated relationships between the response and the explanatory variables.
For example, the estimated difference between Crews C and B can be obtained directly from the original Crew A-reference model:
\[ \begin{aligned} \widehat{\beta}_{C}-\widehat{\beta}_{B} &=0.0144514-(-0.0476194)\\ &=0.0620708. \end{aligned} \]
This is exactly the Crew C coefficient obtained when Crew B is selected as the reference.
The same principle applies to the remaining crew comparisons.
Because the underlying model is unchanged, both parameterizations produce identical:
- Fitted values and residuals
- Residual sum of squares
- Residual standard error
- \(R^2\) and adjusted \(R^2\)
- Overall model F-statistic
- Model-level diagnostic results based on the same fitted values, residuals, and predictor space
Consequently, the diagnostic analyses performed on the original model with Crew A as the reference are equally applicable to the reparameterized model with Crew B as the reference.
This does not imply that the assumptions are necessarily satisfied; their adequacy must still be assessed through the relevant diagnostics.
Model Diagnostics
Visual Assessments
Before conducting formal diagnostic tests, we examined the fitted values and residuals from Model 2 to assess the adequacy of the regression specification.
Since the model contains 16 estimated parameters, including the intercept, it is not possible to visualize the complete relationship between the response and all explanatory variables in a conventional two-dimensional graph.
Instead, we examined the relationship between observed and fitted values, the distribution of residuals against fitted values, and the distribution of residuals relative to a theoretical normal distribution.
The visual assessments presented below apply to both parameterizations of Model 2, since changing the reference crew does not affect fitted values or residuals.
Observed versus Fitted Values
The figure compares the observed log-hours ratio, \(Y_i\), with its fitted value, \(\widehat{Y}_i\), for all 300 projects.
Although the model includes multiple explanatory variables, each observation is represented using only two coordinates:
\[ (\widehat{Y}_i,Y_i). \]
The dashed diagonal line represents perfect agreement between observed and fitted values:
\[ Y_i=\widehat{Y}_i. \]
Observations above this line have positive residuals, indicating that their actual-to-estimated labour ratios were higher than predicted by the model. Conversely, observations below the line have negative residuals.
We observe a positive relationship between fitted and observed log-hours ratios. However, there is considerable dispersion around the diagonal line, suggesting that the included explanatory variables do not account for all project-level variation.
This is consistent with the model’s coefficient of determination, \(R^2=0.3075\), indicating that approximately 30.8% of the variation in the observed log-hours ratio is explained by the fitted regression.
Several observations also lie substantially above the diagonal line, indicating projects with unusually high actual labour requirements relative to the model’s fitted values.
These observations warrant further investigation and are examined later using studentized residuals, leverage, and Cook’s distance.
Residuals versus Fitted Values
The residual for observation \(i\) is defined as
\[ e_i=Y_i-\widehat{Y}_i. \]
In our plot, the residuals appear to be distributed on both sides of zero, with no obvious curvature or consistent increasing or decreasing trend.
Furthermore, we do not observe a pronounced funnel-shaped pattern, which would otherwise suggest that residual variance increases or decreases systematically with fitted values.
These observations provide some visual support for the chosen functional form and do not reveal an obvious violation of the constant-variance assumption.
Formal specification and heteroskedasticity tests are therefore considered in the following sections.
Normal Q-Q Plot
The asymmetric departures, particularly in the upper tail, suggest positive skewness or a heavier-than-normal right tail in the residual distribution.
Consequently, the Q-Q plot raises concerns about the normality assumption underlying classical small-sample regression inference.
We investigate this further through formal residual-normality tests.
Stem-and-Leaf Plot
As an additional visual assessment, we constructed a stem-and-leaf plot of the regression residuals using the following R command:
stem(resid(model_2), scale = 2)The resulting output was:
The decimal point is 2 digit(s) to the left of the |
-18 | 7
-17 | 20
-16 |
-15 | 84
-14 | 7
-13 | 4
-12 | 8760
-11 | 88553
-10 | 882110
-9 | 86432111
-8 | 98664320
-7 | 9987333
-6 | 997776654333200
-5 | 9888665544322210
-4 | 999966543321100
-3 | 9999877776433321100000
-2 | 77755443222110
-1 | 887766543322110
-0 | 986665444331110
0 | 01111223446666788999
1 | 002233334455667779
2 | 002345678
3 | 0013356677889
4 | 00111122345666678899
5 | 012268
6 | 05566899
7 | 0012366789
8 | 0015556
9 | 0224669
10 | 0378
11 | 345
12 | 3
13 | 12489
14 | 0
15 | 39
16 | 8
17 |
18 |
19 | 366
20 | 5
21 |
22 |
23 | 356
24 |
25 | 0This asymmetry, together with the relatively large positive observations in the upper tail, is consistent with the departure from normality observed in the Q-Q plot.
The stem-and-leaf plot therefore provides further descriptive evidence of a more pronounced positive residual tail.
Although these observations suggest that the residual distribution is not perfectly symmetric, the plot alone is insufficient to determine whether the departures from normality are statistically significant.
Formal Statistical Diagnostics
Functional Form: Ramsey RESET Test
To investigate potential functional-form misspecification, we performed the Ramsey Regression Equation Specification Error Test (RESET). (Ramsey 1969)
The hypotheses for the Ramsey RESET test are:
- \(H_0\): The model does not exhibit the functional-form misspecification detected by the additional RESET terms.
- \(H_1\): The model exhibits evidence of functional-form misspecification.
The test was performed in R using:
lmtest::resettest(model_2)Results
| Statistic | Value |
|---|---|
| Test | Ramsey RESET |
| F-statistic | 2.8734 |
| Numerator degrees of freedom | 2 |
| Denominator degrees of freedom | 282 |
| p-value | 0.05816 |
| Significance level | 0.05 |
Since
\[ p=0.05816>0.05, \]
we fail to reject the null hypothesis at the \(\alpha = 5\%\) significance level.
Therefore, the RESET test does not provide sufficient statistical evidence to conclude that the original model is misspecified at this significance level.
However, the p-value is relatively close to 0.05. At a 10% significance level, the null hypothesis would be rejected.
The residuals-versus-fitted plot presented earlier did not reveal an obvious systematic curvature, which is broadly consistent with the absence of strong evidence against the specification at the 5% level.
Residual Independence: Durbin-Watson Test
We examined first-order autocorrelation using the Durbin-Watson test (Durbin and Watson 1950).
The hypotheses are:
- \(H_0\): There is no first-order autocorrelation (\(\rho=0\)).
- \(H_1\): There is positive first-order autocorrelation (\(\rho>0\)).
The test was performed in R using:
lmtest::dwtest(model_2)Results
| Statistic | Value |
|---|---|
| Test | Durbin-Watson |
| DW statistic | 1.9166 |
| p-value | 0.2338 |
| Alternative hypothesis | Positive autocorrelation |
| Significance level | 0.05 |
Since \(p=0.2338>0.05\), we fail to reject the null hypothesis at the 5% significance level.
Furthermore, the Durbin-Watson statistic of 1.9166 is relatively close to 2, suggesting little evidence of positive first-order autocorrelation in the observation sequence used by the test.
Normality: Shapiro-Wilk Test
We assessed residual normality using the Shapiro-Wilk test (Shapiro and Wilk 1965).
The hypotheses are:
- \(H_0\): The residuals follow a normal distribution.
- \(H_1\): The residuals do not follow a normal distribution.
The test was performed in R using:
shapiro.test(residuals(model_2))Results
| Statistic | Value |
|---|---|
| Test | Shapiro-Wilk |
| W-statistic | 0.98042 |
| p-value | 0.0004005 |
| Significance level | 0.05 |
Since
\[ p=0.0004005<0.05, \]
we reject the null hypothesis at the 5% significance level.
However, statistical significance does not necessarily imply that the departure from normality is substantial enough to materially affect regression inference.
With a sample size of 300 observations, the Shapiro-Wilk test can detect relatively small deviations from normality. Our Q-Q plot showed that the central portion of the residual distribution followed the theoretical normal quantiles relatively closely. The principal departures occurred in the upper tail, where several unusually large positive residuals were observed. Similarly, the stem-and-leaf plot suggested a more pronounced positive tail. These observations suggest that the rejection of normality may be influenced by a relatively small number of unusually large positive residuals, rather than a substantial departure throughout the entire distribution.
Heteroskedasticity: Breusch-Pagan Test
We assessed heteroskedasticity using the studentized Breusch-Pagan test (Breusch and Pagan 1979; Koenker 1981).
The hypotheses are:
- \(H_0\): The error terms are homoskedastic (constant conditional variance).
- \(H_1\): The error terms exhibit heteroskedasticity (nonconstant conditional variance).
The test was performed in R using:
lmtest::bptest(model_2)Results
| Statistic | Value |
|---|---|
| Test | Studentized Breusch-Pagan |
| BP statistic | 18.058 |
| Degrees of freedom | 15 |
| p-value | 0.2596 |
| Significance level | 0.05 |
Since
\[ p=0.2596>0.05, \]
we fail to reject the null hypothesis at the 5% significance level.
Therefore, the Breusch-Pagan test does not provide sufficient statistical evidence of heteroskedasticity in Model 2.
This finding is consistent with our earlier visual assessment of the residuals-versus-fitted plot, where the residuals appeared to be scattered around zero without a pronounced funnel-shaped pattern or systematic change in dispersion.
Multicollinearity: Variance Inflation Factor (VIF)
We evaluated multicollinearity using generalized variance inflation factors (Fox and Monette 1992).
To assess multicollinearity among the explanatory variables, we calculated the Variance Inflation Factor (VIF) using:
car::vif(model_2)Results
| Variable | GVIF | Df | Adjusted GVIF |
|---|---|---|---|
| Crew | 1.0863 | 3 | 1.0139 |
| Year | 1.0258 | 1 | 1.0128 |
| Log(Patio area) | 1.0474 | 1 | 1.0234 |
| Skid-steer access | 1.0527 | 1 | 1.0260 |
| Replacement | 1.0455 | 1 | 1.0225 |
| Cutting complexity | 1.0266 | 1 | 1.0132 |
| Against house | 1.0436 | 1 | 1.0215 |
| Restoration method | 1.0475 | 1 | 1.0235 |
| Difficult soil | 1.0304 | 1 | 1.0151 |
| Utility conflict | 1.0264 | 1 | 1.0131 |
| Delivery issue | 1.0387 | 1 | 1.0192 |
| Wet-weather fraction | 1.0474 | 1 | 1.0234 |
| Hot-weather fraction | 1.0363 | 1 | 1.0180 |
Note: Adjusted GVIF = \(\mathrm{GVIF}^{1/(2Df)}\), accounting for the degrees of freedom associated with each predictor.
All GVIF values are close to 1, with the largest being 1.0863 for crew assignment. The adjusted GVIF values are also close to 1.
These results indicate no evidence of substantial multicollinearity among the explanatory variables included in Model 2.
Consequently, multicollinearity does not appear to be a concern for the stability or precision of our regression coefficient estimates.
Heteroskedasticity-Robust Inference
We calculated HC3 heteroskedasticity-robust standard errors (MacKinnon and White 1985), implemented using the R package sandwich (Zeileis 2004).
Following the diagnostic analysis, we calculated HC3 heteroskedasticity-robust standard errors to assess the statistical significance of our regression coefficients.
Although the Breusch-Pagan test did not provide significant evidence of heteroskedasticity, HC3 inference was used as a precaution against potential nonconstant error variance.
HC3 modifies the estimated standard errors without changing the OLS coefficients. It does not correct model misspecification or dependence between observations.
The following results use Crew B as the reference category.
Adjusted Crew Comparisons
| Crew comparison | Coefficient | Adjusted difference | HC3 p-value |
|---|---|---|---|
| A relative to B | 0.047619 | +4.88% | 0.000173 |
| C relative to B | 0.062071 | +6.40% | \(4.02\times10^{-7}\) |
| D relative to B | 0.100957 | +10.62% | \(1.64\times10^{-13}\) |
All three comparisons were statistically significant at the 5% level under HC3 inference.
These findings indicate that Crews A, C, and D had higher actual-to-estimated labour ratios than Crew B, after accounting for the measured explanatory variables.
Operational Predictors
| Predictor | Coefficient | Adjusted difference | HC3 p-value |
|---|---|---|---|
| Difficult soil | 0.043264 | +4.42% | 0.000136 |
| Delivery issue | 0.052064 | +5.34% | 0.01334 |
| Wet-weather fraction | 0.093276 | +0.94% per 0.10 increase | \(8.80\times10^{-6}\) |
| Utility conflict | 0.039300 | +4.01% | 0.05179 |
Difficult soil, delivery issues, and wet-weather exposure were significantly associated with higher actual-to-estimated labour ratios.
However, utility conflicts were significant under conventional OLS inference (\(p=0.0206\)) but not under HC3 inference (\(p=0.0518\)).
Consequently, the evidence for an association between utility conflicts and labour overruns is less conclusive.
These results are consistent with the findings presented in the companion management decision brief. However, the coefficients represent statistical associations rather than causal effects.
Influential Observations
We used Cook’s distance to identify potentially influential observations (Cook 1977).
We examined Model 2 for unusual observations using externally studentized residuals, leverage values, and Cook’s distance.
These diagnostics identify different characteristics of an observation. A project may have an unusually large residual without having high leverage, while a high-leverage project may be fitted accurately and have little influence on the regression estimates.
Outliers in the Response Variable
We first identified projects with unusually large residuals using externally studentized residuals.
An observation was flagged when:
\[ |t_i|>3, \]
where \(t_i\) denotes the externally studentized residual for project \(i\).
Four projects exceeded this threshold:
| Project | Crew | Estimated hours | Actual hours | Original overrun | Studentized residual |
|---|---|---|---|---|---|
| P050 (2025) | C | 9.55 | 12.60 | +31.6% | 3.48 |
| P194 (2024) | A | 12.30 | 17.10 | +38.3% | 3.23 |
| P076 (2025) | D | 14.30 | 19.10 | +33.8% | 3.19 |
| P029 (2025) | B | 9.32 | 11.20 | +20.1% | 3.13 |
Note: Labour hours are rounded for presentation; overruns were calculated using unrounded values.
All four observations had positive studentized residuals, indicating that their actual-to-estimated labour ratios were substantially higher than those predicted by Model 2, even after accounting for their recorded project characteristics and operational conditions.
Interestingly, three of these projects (P050, P076, and P029) had no recorded difficult soil, utility conflicts, or delivery issues. P194 had a recorded utility conflict.
These findings suggest that certain labour overruns may be associated with unrecorded complications or other factors not captured by the model. However, the residuals alone cannot establish their causes.
High-Leverage Observations
We next examined leverage, which measures how unusual an observation’s combination of explanatory variables is relative to the rest of the dataset.
For Model 2, we have \(n=300\) observations and \(p=16\) estimated coefficients, including the intercept.
We used the screening threshold:
\[ h_{ii}>\frac{2p}{n} =\frac{2(16)}{300} \approx0.1067. \]
Four projects exceeded this threshold:
| Project | Crew | Patio area (sq ft) | Leverage | Cook’s distance |
|---|---|---|---|---|
| P052 (2025) | B | 760 | 0.133 | 0.000019 |
| P046 (2025) | C | 232 | 0.124 | 0.00227 |
| P163 (2024) | D | 453 | 0.119 | 0.0381 |
| P103 (2025) | C | 426 | 0.111 | 0.00800 |
All four projects had recorded utility conflicts, indicating that they shared a relatively uncommon combination of explanatory variables. None exceeded the more restrictive leverage threshold of \(3p/n=0.16\).
Importantly, high leverage does not necessarily imply an unusually large residual or substantial influence on the regression coefficients.
For example, P052 had the highest leverage but was fitted almost exactly by the model, resulting in a Cook’s distance of approximately 0.000019.
In contrast, P163 combined high leverage with a relatively large residual and a Cook’s distance of 0.0381. This project experienced both a utility conflict and a delivery issue, requiring approximately 41.5% more labour than originally estimated.
Consequently, P163 warrants closer investigation, despite not meeting the \(|t_i|>3\) response-outlier threshold.
Visual Identification of Unusual Observations
The figure highlights the four response outliers in red and the four high-leverage observations in green.
The distinction is important. Response outliers are located relatively far above the diagonal fitted-value line, whereas high-leverage observations are not necessarily associated with large prediction discrepancies.
P052 illustrates this distinction particularly well: despite its high leverage, it lies almost directly on the fitted-value line.
P163, however, lies noticeably above the line and is also a high-leverage observation, making it a stronger candidate for investigating potential influence on the regression estimates.
Influential Observations: Cook’s Distance
Finally, we used Cook’s distance to identify observations that could substantially affect the estimated regression coefficients.
An observation was flagged for further investigation when:
\[ D_i>\frac{4}{n} =\frac{4}{300} \approx0.0133. \]
Of the 300 observations, 18 exceeded this threshold.
The ten projects with the largest Cook’s distances were:
| Project | Crew | Cook’s distance | Studentized residual |
|---|---|---|---|
| P050 | C | 0.0710 | +3.48 |
| P194 | A | 0.0500 | +3.23 |
| P207 | D | 0.0403 | +2.70 |
| P138 | D | 0.0394 | +2.66 |
| P163 | D | 0.0381 | +2.14 |
| P076 | D | 0.0256 | +3.19 |
| P174 | B | 0.0228 | -2.11 |
| P110 | D | 0.0223 | -2.53 |
| P084 | B | 0.0210 | -2.14 |
| P140 | D | 0.0203 | +1.89 |
Notably, six of the ten observations with the highest Cook’s distances belonged to Crew D.
This is particularly relevant because Crew D had the largest adjusted labour ratio relative to Crew B in our regression analysis.
These observations raise an important question: Is Crew D’s estimated coefficient being disproportionately influenced by a relatively small number of projects?
We investigate this possibility in the subsequent sensitivity analysis by refitting the model after excluding selected Crew D observations.
Cook’s distance identifies observations that warrant investigation; exceeding the threshold does not establish that a project contains an error or should be removed.
Summary
The influence diagnostics provide the statistical basis for the five projects recommended for management review:
- P050, P194, P076, and P029: Unusually large positive studentized residuals.
- P163: High leverage, substantial original labour overrun, multiple recorded complications, and a Cook’s-distance flag.
All five projects also exceeded the Cook’s-distance screening threshold.
These observations should therefore be investigated, but they remain in the primary regression model unless a valid reason for exclusion is established.
Sensitivity Analysis
The Cook’s distance analysis identified 18 potentially influential observations, including six Crew D projects among the ten highest Cook’s distances.
Since Crew D exhibited the largest adjusted labour overrun relative to Crew B, we conducted a sensitivity analysis to determine whether these observations disproportionately influenced its estimated regression coefficient.
Individual Exclusion of Selected Projects
We first refitted Model 2 six times, excluding one selected Crew D project at a time.
The original regression used Crew A as the reference category, with an estimated Crew D coefficient of \(\widehat{\beta}_D=0.053337\).
| Project removed | Cook’s distance | Refitted Crew D coefficient | Change from original |
|---|---|---|---|
| P207 | 0.0403 | 0.049920 | -0.003417 |
| P138 | 0.0394 | 0.051045 | -0.002292 |
| P163 | 0.0381 | 0.051746 | -0.001591 |
| P076 | 0.0256 | 0.050511 | -0.002826 |
| P110 | 0.0223 | 0.055306 | +0.001969 |
| P140 | 0.0203 | 0.051720 | -0.001617 |
The coefficient remained positive in every refitted model, ranging from approximately 0.0499 to 0.0553.
Removing five of the six projects decreased the coefficient, while removing P110 increased it. This is consistent with P110 having a negative studentized residual.
These results suggest that no single observation among the six selected projects was responsible for Crew D’s positive adjusted coefficient.
Simultaneous Exclusion of Six Projects
We subsequently removed all six selected Crew D projects simultaneously and refitted Model 2 using the remaining 294 observations.
Regression coefficients with Crew A as the reference
| Coefficient | Full model (n=300) | Without six (n=294) |
|---|---|---|
| Crew B | -0.047619 | -0.046647 |
| Crew C | +0.014451 | +0.014928 |
| Crew D | +0.053337 | +0.043050 |
Crew D’s coefficient decreased from 0.053337 to 0.043050, while the coefficients for Crews B and C remained relatively stable.
Converting these coefficients into percentage differences using
\[ 100(e^{\widehat{\beta}}-1), \]
the adjusted difference between Crews D and A decreased from 5.48% to 4.40%.
Nevertheless, Crew D’s coefficient remained statistically significant under HC3 robust inference (\(p=0.0009854\)).
Sensitivity Analysis Using Crew B as the Reference
Since Crew B was selected as the reference category in the management report, we also expressed the sensitivity results relative to Crew B.
Adjusted crew comparisons
| Crew comparison | Full model (n=300) | Without six (n=294) | HC3 p-value without six |
|---|---|---|---|
| A relative to B | +4.88% | +4.78% | 0.000194 |
| C relative to B | +6.40% | +6.35% | \(5.20\times10^{-7}\) |
| D relative to B | +10.62% | +9.38% | \(1.52\times10^{-13}\) |
The adjusted difference between Crews D and B decreased from 10.62% to 9.38%, a reduction of approximately 1.24 percentage points.
However, the difference remained statistically significant under HC3 robust inference.
The adjusted differences for Crews A and C remained relatively stable.
Therefore, although the six selected influential observations affected the magnitude of Crew D’s estimated labour discrepancy, they did not account for the entire difference relative to Crew B.
Sensitivity of Operational Predictors
We also examined whether excluding these projects affected the coefficients associated with operational complications.
| Predictor | Full model coefficient | Without six | HC3 p-value without six |
|---|---|---|---|
| Difficult soil | 0.043264 | 0.045667 | 0.000021 |
| Utility conflict | 0.039300 | 0.025610 | 0.164847 |
| Delivery issue | 0.052064 | 0.035422 | 0.050104 |
| Wet-weather fraction | 0.093276 | 0.091394 | 0.00000665 |
The difficult-soil and wet-weather coefficients remained relatively stable.
However, the utility-conflict and delivery-issue coefficients decreased noticeably. The associated HC3 p-values also increased, suggesting that the statistical evidence for these relationships is more sensitive to the selected observations.
Conclusion
The sensitivity analysis indicates that Crew D’s adjusted labour discrepancy is not driven entirely by the six selected influential observations, although these projects contribute to its estimated magnitude.
Importantly, the adjusted comparisons relative to Crew B remained positive and statistically significant after the exclusions, supporting the stability of the principal crew-level findings presented in the management report.
Nevertheless, influential observations are not necessarily erroneous. Since there is no established reason to exclude these projects, all 300 observations were retained in the primary regression model.
The reduced-sample analysis is presented strictly as a sensitivity check rather than an alternative preferred specification.
Conclusion: Model Reliability and Intended Use
The objective of Model 2 was primarily explanatory rather than predictive. We sought to investigate which measured project characteristics and operational conditions were associated with labour estimating discrepancies, rather than to predict the precise labour requirements of future installations.
The model achieved an \(R^2\) of 0.3075, explaining approximately 30.8% of the observed variation in the log actual-to-estimated labour ratio. Although this indicates that substantial project-level variation remains unexplained, it does not necessarily undermine the model’s usefulness for identifying statistical associations.
Our diagnostic analysis produced the following findings:
- Functional form: The Ramsey RESET test did not reject the specification at the 5% level (\(p=0.0582\)), although some uncertainty remains.
- Homoskedasticity: The Breusch-Pagan test found no statistically significant evidence of heteroskedasticity (\(p=0.2596\)).
- Normality: The Shapiro-Wilk test rejected residual normality (\(p=0.0004\)). Visual diagnostics suggested that departures were concentrated primarily in the upper tail.
- Multicollinearity: Variance inflation factors were close to one, indicating no substantial multicollinearity among the included predictors.
- Independence: The Durbin-Watson test did not identify significant positive first-order autocorrelation in the supplied observation order. However, independence between projects, particularly those completed by the same crew, was not established.
The influence diagnostics identified several unusual observations. Nevertheless, the sensitivity analysis demonstrated that the principal crew-level findings were not entirely dependent on the six selected influential Crew D projects. In particular, Crew D’s adjusted labour ratio remained significantly higher than Crew B’s after their exclusion.
Taken together, these findings suggest that the model provides useful, but qualified, evidence of associations between labour estimating discrepancies, crew assignment, and recorded operational conditions.
However, the model has limitations. Important characteristics such as foreman experience, workforce composition, and the severity of operational complications were not recorded. Consequently, the estimated regression coefficients cannot be interpreted as causal effects.
Furthermore, we did not assess the model’s predictive performance on unseen projects using out-of-sample validation. Its predictive reliability therefore remains unknown, and the relatively modest in-sample \(R^2\) does not justify using it to forecast individual project labour requirements without further evaluation.
Overall, Model 2 is best suited as an explanatory tool to support management investigations. It identifies recurring associations and unusual projects that warrant attention, but it should not be used to assign blame to individual crews or automatically revise labour estimates.
Finally, because the project records are synthetic, the findings demonstrate the application of the statistical methodology rather than establishing conclusions about the performance of an actual landscaping contractor.