Set-up and background

The assignment is worth 100 points. There are 13 questions. You should have the following packages installed:

library(tidyverse)
library(patchwork)
library(fixest)

In this problem set you will replicate some results from the paper “Do Workers Work More if Wages Are High? Evidence from a Randomized Field Experiment” (Fehr and Goette, AER 2007).

1 Big picture

[Q1] Recall the taxi cab studies where reference dependence is studied using observational data. What can an experimental study do that an observational study can’t?

Experimental studies offer insights that observational studies cannot, as researchers are allowed to manipulate variables directly. For example, in the taxi cab studies, an experimental approach might involve setting different income targets for drivers and observing how these controlled changes impact their labor supply decisions. This allows researchers to clarify whether drivers’ behavior is genuinely driven by reference-dependent preferences, rather than simply correlated with them. Furthermore, experimental studies enable the control of confounding variables, such as demand or weather. Observational studie effectively capture natural behavior, while experiments are better suited to test hypotheses about causal relationships, particularly in exploring mechanisms like loss aversion and reference dependence.

[Q2] Summarize the field experiment design.

Fehr and Goette (2007) conducted a field experiment with Veloblitz bike messengers to study the effect of a temporary wage increase on labor supply. The experiment involved a 25% increase in the commision rate for a randomly assigned treatment group, while a control group continued to receive thei regular commission rate. To ensure balanced results, the experimented was split into two periods. In the first period, Group A received the wage increase, while Group B acted as a control.In the second period, the roles reversed with Group B receiving the increase while Group A acting as the control. The cross-over design allowed each group to serve as both treatment and control, canceling out individual differences which ensured that the estimated effect was independent of individual heterogeneity.

The researchers ensured that the wage increase was fully anticipated and temporary with additional earnings disbursed after the study concluded. This delayed payout controlled for income effects enabling a cleaner isolation of the substitution effect. Time effects due to demand variations were managed by monitoring a second messenger service operating in the same market. Data included the number of shifts and the number of deliveries or revenues per shift, allowing the study to assess both hours worked and exerted effort. Thus, the experiment provided a rigorous design to explore how temporary wage increases influence both the labor hours and intensity of work.

[Q3] Summarize the main results of the field experiment. Distinguish between hours and effort.

The field experiment by Fehr and Goette (2007) found significant effects of a temporary wage increase on both hours worked and effort ecerted among Veloblitz bike messengers. The temporary 25% wage increase resulted in a positive effect on the total revenue per messenger, indicating intertemporal elasticity of substitution between 1.12% and 1.25%, which is higher than prior studies reported.

For hours worked, the experiment showed that the messengers worked an average of 4 additional shifts during the treatement period. Wih an average of 11.925 shifts during the treatement period, the wage elasticity of shifts was estimated between 11.2% and 14.1%. This responsiveness to the wage increase in shifts was greater than the response in total reveue, suggesting a notable increase in the number of shifts worked when wages rose.

In contrast, effort per shift declined by roughly 6% during the wage increase, with a wage elasticity of -24%. This reduction in per-shift revenue suggests that while messengers increased their hours, they may have reduced effort per shift when paid a higher commission. This pattern aligns with a target income model under loss aversion: messengers set daily income targets and adjusted their effort downward once those targets were achievable with fewer deliveries due to the higher commission rate.

Finally, the study observed that this decrease in effort per shift was more pronounced among messengers with higher levels of loss aversion, evidenced by their behavior in lottery tasks, whereas those with lower loss aversion did not significantly reduce efffort. These findings support the idea that loss aversion and reference-dependent preferences play a role in labor supply decisions, as messengers sough to meet specific income targets rather than simply maximizing effort.

2 Replication

Use theme_classic() for all plots.

2.1 Correlations in revenues across firms

For this section please use dailycorrs.csv.

dailycorrs <- read.csv("dailycorrs.csv")

[Q4] The authors show that earnings at Veloblitz and Flash are correlated. Show this with a scatter plot with a regression line and no confidence interval. Title your axes and the plot appropriately. Do not print the plot but assign it to an object called p1.

library(ggplot2)

# Create the scatter plot with regression line and classic theme
p1 <- ggplot(dailycorrs, aes(x = logf, y = logv)) +
  geom_point() +  # Add scatter plot points
  geom_smooth(method = "lm", se = FALSE, color = "blue") +  # Add regression line without confidence interval
  labs(title = "Correlation between Earnings at Veloblitz (v) and Flash (f)",
       x = "Log Earnings at Flash (logf)",
       y = "Log Earnings at Veloblitz (logv)") +
  theme_classic()  # Apply classic theme

p1

[Q5] Next plot the kernel density estimates of revenues for both companies. Overlay the distributions and make the densities transparent so they are easily seen. Title your axes and the plot appropriately. Do not print the plot but assign it to an object called p2.

# Create the kernel density plot for revenues at both companies
p2 <- ggplot(dailycorrs) +
  geom_density(aes(x = logf, fill = "Flash (logf)"), alpha = 0.5) +  # Density for Flash with transparency
  geom_density(aes(x = logv, fill = "Veloblitz (logv)"), alpha = 0.5) +  # Density for Veloblitz with transparency
  labs(title = "Kernel Density Estimates of Revenues for Veloblitz and Flash",
       x = "Log Earnings",
       y = "Density") +
  scale_fill_manual(values = c("Flash (logf)" = "blue", "Veloblitz (logv)" = "red")) +  # Colors for the fills
  theme_classic()  # Apply classic theme

p2

[Q6] Now combine both plots using library(patchwork) and label the plots with letters.

# Combine both plots with labels
combined_plot <- p1 + p2 + 
  plot_annotation(
    title = "Analysis of Earnings Correlation and Revenue Distributions",
    tag_levels = 'a'  # Label plots with letters (a, b)
  )

combined_plot

2.2 Tables 2 and 3

For this section please use tables1to4.csv.

tables1to4 <- read_csv('tables1to4.csv')

2.2.1 Table 2

On page 307 the authors write:

“Table 2 controls for individual fixed effects by showing how, on average, the messengers’ revenues deviate from their person-specific mean revenues. Thus, a positive number here indicates a positive deviation from the person-specific mean; a negative number indicates a negative deviation.”

[Q7] Fixed effects are a way to control for heterogeneity across individuals that is time invariant. Why would we want to control for fixed effects? Give a reason how bike messengers could be different from each other, and how these differences might not vary over time.

Controlling for fixed effects allows researchers to account for individual differences that remain constant over time but might influence the dependent variables. By controlling for these constant characteristics, we can isolate the effect on the primary outcome, in this case, labor supply and effort. This ensures that any changes observed are due to the treatement itself, rather than confounding variables that vary across indivudlas but not over time.

In the case of bike messengers, they may have inherent time invariant characteristics that make them different from each other, such as work habits, motivation levels, personality and physical stamina. These differences do not change over time and could affect their labor supply and productivity. For example, some messengers might be naturally more motivated to work long hours or deliver more packages due to personal goals, while others may prefer shorter shifts or have a lower capacity for effort. Without controlling for these differences, any observed differences in labor supply or effort could be confounded by these underlying characteristics rather than reflecting the true effect of the wage increase.

[Q8] Create a variable called totrev_fe and add it to the dataframe. This requires you to “average out” each individual’s revenue for a block from their average revenue: \(x_i^{fe} = x_{it} - \bar{x}_i\) where \(x_i^{fe}\) is the fixed effect revenue for \(i\).

library(dplyr)

tables1to4 <- tables1to4 %>%
  group_by(fahrer) %>%  # Group by individual 
  mutate(
    avg_rev = mean(totrev, na.rm = TRUE),  # Calculate the average revenue for each individual
    totrev_fe = totrev - avg_rev  # Subtract individual average from the total revenue
  ) %>%
  ungroup()  # Remove grouping after calculation

head(tables1to4)
## # A tibble: 6 × 17
##   vebli block totrev shifts  high firstblock block1 block2 block3 maxhigh other
##   <dbl> <dbl>  <dbl>  <dbl> <dbl>      <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <dbl>
## 1     1     1  3091.     13     0      14836      1      0      0       1     0
## 2     1     2  1210       5     1      14864      0      1      0       1     0
## 3     1     3   782       3     0      14913      0      0      1       1     0
## 4     1     1  2308.      8     0      14836      1      0      0       1     0
## 5     1     2  2148       6     1      14864      0      1      0       1     0
## 6     1     3  1936.      6     0      14913      0      0      1       1     0
## # ℹ 6 more variables: group <chr>, experiment <dbl>, odd <dbl>, fahrer <dbl>,
## #   avg_rev <dbl>, totrev_fe <dbl>

[Q9] Use summarise() to recreate the findings in Table 2 for “Participating Messengers” using your new variable totrev_fe. (You do not have to calculate the differences in means.)

In addition to calculating the fixed-effect controled means, calculate the standard errors. Recall the standard error is \(\frac{s_{jt}}{\sqrt{n_{jt}}}\) where \(s_{jt}\) is the standard deviation for treatment \(j\) in block \(t\) and \(n_{jt}\) are the corresponding number of observations.

(Hint: use n() to count observations.) Each calculation should be named to a new variable. Assign the resulting dataframe to a new dataframe called df_avg_revenue.

# Create a new variable for Group A and Group B based on odd/even messenger ID
df_avg_revenue <- tables1to4 %>%
  filter(maxhigh == 1) %>%  # Only include participating messengers
  mutate(group_type = ifelse(odd %% 2 == 1, "A", "B")) %>%  # Group A for odd, B for even
  group_by(group_type, block) %>%  # Group by both group_type (A/B) and block (1, 2, 3)
  summarise(
    mean_totrev_fe = mean(totrev_fe, na.rm = TRUE),  # Mean of totrev_fe for each block and group
    sd_totrev_fe = sd(totrev_fe, na.rm = TRUE),  # Standard deviation of totrev_fe for each block and group
    n_totrev_fe = n(),  # Number of observations in each block and group
    se_totrev_fe = sd_totrev_fe / sqrt(n_totrev_fe),  # Standard error for each block and group
    .groups = 'drop'  # Remove the grouping after summarising
  )

# View the resulting dataframe
df_avg_revenue
## # A tibble: 6 × 6
##   group_type block mean_totrev_fe sd_totrev_fe n_totrev_fe se_totrev_fe
##   <chr>      <dbl>          <dbl>        <dbl>       <int>        <dbl>
## 1 A              1          -48.9        1654.          21         361.
## 2 A              2          722.          905.          22         193.
## 3 A              3         -675.         1354.          22         289.
## 4 B              1         -120.         1319.          19         303.
## 5 B              2         -278.         1076.          20         241.
## 6 B              3          392.         1121.          20         251.

[Q10] Plot df_avg_revenue. Use points for the means and error bars for standard errors of the means.

To dodge the points and size them appropriately, use

geom_point(position=position_dodge(width=0.5), size=4)

To place the error bars use

geom_errorbar(aes(
  x=block, 
  ymin = [MEAN] - [SE], ymax = [MEAN] + [SE]),
  width = .1,
  position=position_dodge(width=0.5))

You will need to replace [MEAN] with whatever you named your average revenues and [SE] with whatever you named your standard errors.

# Plot the df_avg_revenue
ggplot(df_avg_revenue, aes(x = block, y = mean_totrev_fe, color = group_type, shape = group_type)) +
  geom_point(position = position_dodge(width = 0.5), size = 4) +  # Points for the means
  geom_errorbar(aes(ymin = mean_totrev_fe - se_totrev_fe, ymax = mean_totrev_fe + se_totrev_fe), 
                position = position_dodge(width = 0.5), width = 0.25) +  # Error bars for SE
  labs(x = "Block", y = "Fixed-effect Controlled Mean Revenue (totrev_fe)", 
       title = "Fixed-effect Controlled Means with Standard Errors") + 
  theme_classic() +  # Clean theme
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

2.2.2 Table 3

[Q11] Recreate the point estimates in Model (1) in Table 3 by hand (you don’t need to worry about the standard errors). Assign it to object m1. Recreating this model requires you to control for individual fixed effects and estimate the following equation where \(\text{H}\) is the variable high, \(\text{B2}\) is the second block (block == 2) and \(\text{B3}\) is the third block (block == 3):

\[ y_{ijt} - \bar{y}_{ij} = \beta_1 (\text{H}_{ijt} - \bar{\text{H}}_{ij}) + \beta_2 (\text{B2}_{ijt} - \bar{\text{B2}}_{ij}) + \beta_3 (\text{B3}_{ijt} - \bar{\text{B3}}_{ij}) + (\varepsilon_{ijt} - \bar{\varepsilon}_{ij}) \]

df_participating <- tables1to4 %>%
  filter(group %in% c("CG Veloblitz", "TG Veloblitz")) %>%
  group_by(fahrer) %>%
  mutate(
    mean_totrev = mean(totrev),   # Individual mean for totrev
    mean_high = mean(high),       # Individual mean for high
    mean_block2 = mean(block2),   # Individual mean for block2
    mean_block3 = mean(block3)    # Individual mean for block3
  ) %>%
  ungroup() 

df_participating <- df_participating %>%
  mutate(
    y_dif = totrev - mean_totrev,     # Difference of totrev from individual mean
    H_dif = high - mean_high,         # Difference of high from individual mean
    B2_dif = block2 - mean_block2,    # Difference of block2 from individual mean
    B3_dif = block3 - mean_block3     # Difference of block3 from individual mean
  )

# Calculate covariances of each independent variable with y_dif
cov_H_y <- cov(df_participating$H_dif, df_participating$y_dif)
cov_B2_y <- cov(df_participating$B2_dif, df_participating$y_dif)
cov_B3_y <- cov(df_participating$B3_dif, df_participating$y_dif)

# Calculate variances of each independent variable
var_H <- var(df_participating$H_dif)
var_B2 <- var(df_participating$B2_dif)
var_B3 <- var(df_participating$B3_dif)

# Calculate the regression coefficients
beta_H <- cov_H_y / var_H
beta_B2 <- cov_B2_y / var_B2
beta_B3 <- cov_B3_y / var_B3

# Combine results into a named vector for easy viewing
m1 <- c(beta_H = beta_H, beta_B2 = beta_B2, beta_B3 = beta_B3)
m1
##    beta_H   beta_B2   beta_B3 
##  857.3855  373.1783 -253.7205

[Q12] Now recreate the same point estimates using lm and assign it to object m2. You are estimating the model below where where \(\text{F}_i\) is the dummy variable for each messenger (fahrer). Make sure to cluster the standard errors at the messenger level. (Use lmtest and sandwhich for this.)

\[ y_{ijt} = \beta_0 + \beta_1 \text{H}_{ijt} + \beta_2 \text{B2}_{ijt} + \beta_3 \text{B3}_{ijt} + \sum_{i=1}^{n} \alpha_i \text{F}_i + \varepsilon_{ijt} \]

library(lmtest)
library(sandwich)

# Fit the model with individual fixed effects by including fahrer as a factor
m2 <- lm(
  totrev ~ high + block2 + block3 + factor(fahrer), 
  data = df_participating
)

# Clustered standard errors at the messenger level
clustered_se <- vcovCL(m2, cluster = ~ fahrer)

# Output model summary with clustered standard errors
m2_summary <- coeftest(m2, vcov. = clustered_se)
m2_summary
## 
## t test of coefficients:
## 
##                     Estimate  Std. Error     t value  Pr(>|t|)    
## (Intercept)       1.6117e+03  2.9297e+02  5.5012e+00 4.509e-07 ***
## high              1.0336e+03  3.2685e+02  3.1621e+00  0.002222 ** 
## block2           -2.1097e+02  4.9725e+02 -4.2430e-01  0.672516    
## block3           -5.7471e+02  5.4568e+02 -1.0532e+00  0.295454    
## factor(fahrer)2   4.3670e+02  6.0834e-12  7.1785e+13 < 2.2e-16 ***
## factor(fahrer)3   7.4720e+02  6.1385e-12  1.2172e+14 < 2.2e-16 ***
## factor(fahrer)5  -1.2704e+03  1.4648e+02 -8.6727e+00 4.240e-13 ***
## factor(fahrer)6   1.5375e+03  6.1857e-12  2.4856e+14 < 2.2e-16 ***
## factor(fahrer)7   7.5017e+02  7.1903e-12  1.0433e+14 < 2.2e-16 ***
## factor(fahrer)8  -1.0697e+02  6.3490e-12 -1.6848e+13 < 2.2e-16 ***
## factor(fahrer)9   3.9122e+03  6.1472e-12  6.3642e+14 < 2.2e-16 ***
## factor(fahrer)14  8.8400e+02  6.8722e-12  1.2863e+14 < 2.2e-16 ***
## factor(fahrer)15 -3.9967e+01  6.1933e-12 -6.4532e+12 < 2.2e-16 ***
## factor(fahrer)18  1.3274e+03  6.2285e-12  2.1311e+14 < 2.2e-16 ***
## factor(fahrer)19  2.6590e+03  8.7866e-12  3.0262e+14 < 2.2e-16 ***
## factor(fahrer)21 -5.1413e+02  6.1467e-12 -8.3644e+13 < 2.2e-16 ***
## factor(fahrer)22  2.9556e+03  6.3208e-12  4.6760e+14 < 2.2e-16 ***
## factor(fahrer)23  3.3302e+03  7.1260e-12  4.6733e+14 < 2.2e-16 ***
## factor(fahrer)24  1.7592e+03  6.2055e-12  2.8348e+14 < 2.2e-16 ***
## factor(fahrer)25  6.1167e+01  6.2755e-12  9.7468e+12 < 2.2e-16 ***
## factor(fahrer)28 -1.5017e+02  6.2625e-12 -2.3979e+13 < 2.2e-16 ***
## factor(fahrer)30  1.5575e+03  6.1501e-12  2.5325e+14 < 2.2e-16 ***
## factor(fahrer)31 -7.3797e+02  6.2464e-12 -1.1814e+14 < 2.2e-16 ***
## factor(fahrer)32 -5.2610e+02  6.1694e-12 -8.5276e+13 < 2.2e-16 ***
## factor(fahrer)33  2.1732e+03  6.1551e-12  3.5307e+14 < 2.2e-16 ***
## factor(fahrer)34  9.1070e+03  7.1333e-12  1.2767e+15 < 2.2e-16 ***
## factor(fahrer)35 -6.4000e+01  6.1574e-12 -1.0394e+13 < 2.2e-16 ***
## factor(fahrer)36  2.0549e+03  6.1560e-12  3.3380e+14 < 2.2e-16 ***
## factor(fahrer)37  2.1940e+03  6.1419e-12  3.5721e+14 < 2.2e-16 ***
## factor(fahrer)38  1.8369e+03  6.3206e-12  2.9063e+14 < 2.2e-16 ***
## factor(fahrer)42  3.9380e+03  6.1411e-12  6.4126e+14 < 2.2e-16 ***
## factor(fahrer)44  5.2187e+02  6.1318e-12  8.5108e+13 < 2.2e-16 ***
## factor(fahrer)45  5.6977e+02  6.4151e-12  8.8816e+13 < 2.2e-16 ***
## factor(fahrer)49  1.0442e+03  6.2792e-12  1.6629e+14 < 2.2e-16 ***
## factor(fahrer)50  6.5907e+03  6.1428e-12  1.0729e+15 < 2.2e-16 ***
## factor(fahrer)51  5.2890e+03  6.1418e-12  8.6114e+14 < 2.2e-16 ***
## factor(fahrer)52  4.2783e+03  6.1322e-12  6.9768e+14 < 2.2e-16 ***
## factor(fahrer)53  3.4330e+03  6.1367e-12  5.5942e+14 < 2.2e-16 ***
## factor(fahrer)55  9.7627e+02  6.2139e-12  1.5711e+14 < 2.2e-16 ***
## factor(fahrer)56  1.7618e+03  6.2240e-12  2.8307e+14 < 2.2e-16 ***
## factor(fahrer)57  2.4715e+03  6.1432e-12  4.0232e+14 < 2.2e-16 ***
## factor(fahrer)58  1.3532e+03  6.3263e-12  2.1390e+14 < 2.2e-16 ***
## factor(fahrer)60  9.3103e+02  6.1323e-12  1.5182e+14 < 2.2e-16 ***
## factor(fahrer)61  6.6470e+02  6.1854e-12  1.0746e+14 < 2.2e-16 ***
## factor(fahrer)63 -2.6111e+02  1.4648e+02 -1.7825e+00  0.078504 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

[Q13] Now use feols to recreate Model (1), including the standard errors. Assign your estimates to the object m3. You are estimating the model below where where \(\alpha_i\) is the individual intercept (i.e. the individual fixed effect):

\[ y_{ijt} = \alpha_i + \beta_1 \text{H}_{ijt} + \beta_2 \text{B2}_{ijt} + \beta_3 \text{B3}_{ijt} + \varepsilon_{ijt} \]

library(fixest)

# Fit the model with clustered standard errors at the 'fahrer' level
m3 <- feols(
  totrev ~ high + block2 + block3 | fahrer,
  data = df_participating,
  cluster = "fahrer" 
)

summary(m3)
## OLS estimation, Dep. Var.: totrev
## Observations: 124
## Fixed-effects: fahrer: 42
## Standard-errors: Clustered (fahrer) 
##        Estimate Std. Error   t value   Pr(>|t|)    
## high   1033.560    265.202  3.897256 0.00035253 ***
## block2 -210.973    403.458 -0.522911 0.60385066    
## block3 -574.713    442.748 -1.298057 0.20152287    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1,230.2     Adj. R2: 0.59584 
##                 Within R2: 0.123224