Case study: windmill

The windmill data

• Engineer: does amount of electricity generated by windmill depend on how strongly wind blowing?
• Measurements of wind speed and DC current generated at various times.
• Assume the “various times” to be randomly selected — aim to generalize to “this windmill at all times”.
• Research questions:
  • Relationship between wind speed and current generated?
  • If so, what kind of relationship?
  • Can we model relationship to do predictions?

Packages for this section

library(tidyverse)
library(broom)

Reading in the data

my_url <- 
  "http://ritsokiguess.site/datafiles/windmill.csv"
windmill <- read_csv(my_url)
windmill

Strategy

• Two quantitative variables, looking for relationship: regression methods.
• Start with picture (scatterplot).
• Fit models and do model checking, fixing up things as necessary.
• Scatterplot:
  • 2 variables, DC_output and wind_velocity.
  • First is output/response, other is input/explanatory.
  • Put DC_output on vertical scale.
• Add trend, but don’t want to assume linear:

ggplot(windmill, aes(y = DC_output, x = wind_velocity)) +
  geom_point() + geom_smooth(se = FALSE) 

Scatterplot

Comments

• Definitely a relationship: as wind velocity increases, so does DC output. (As you’d expect.)
• Is relationship linear? To help judge, geom_smooth smooths scatterplot trend.
• Trend more or less linear for while, then curves downwards (levelling off?).
• Straight line probably not so good here.

Fit a straight line (and see what happens)

DC.1 <- lm(DC_output ~ wind_velocity, data = windmill)
summary(DC.1)

Call:
lm(formula = DC_output ~ wind_velocity, data = windmill)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.59869 -0.14099  0.06059  0.17262  0.32184 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.13088    0.12599   1.039     0.31    
wind_velocity  0.24115    0.01905  12.659 7.55e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2361 on 23 degrees of freedom
Multiple R-squared:  0.8745,    Adjusted R-squared:  0.869 
F-statistic: 160.3 on 1 and 23 DF,  p-value: 7.546e-12

Another way of looking at the output

• The standard output tends to go off the bottom of the page rather easily. Package broom has these:

glance(DC.1)

showing R-squared, and:

tidy(DC.1)

showing intercept and slope and their significance.

Comments

• Strategy: lm actually fits the regression. Store results in a variable. Then look at the results, eg. via summary or glance/tidy.
• My strategy for model names: base on response variable (or data frame name) and a number. Allows me to fit several models to same data and keep track of which is which.
• The outputs from glance and tidy are dataframes, so are helpful for pulling out results from, eg.:

glance(DC.1) %>% 
  pluck("r.squared", 1) %>% 
  round(3) -> rsq1

• Results actually pretty good: wind.velocity strongly significant, R-squared (0.874) high.
• How to check whether regression is appropriate? Look at the residuals, observed minus predicted, plotted against fitted (predicted).

Scatterplot, but with line

ggplot(windmill, aes(y = DC_output, x = wind_velocity)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)

… alternatively

ggplot(windmill, aes(y = DC_output, x = wind_velocity)) +
  geom_point() + geom_smooth(method = "lm")

To plot residuals

• First make dataframe containing regression information (residuals, fitted values etc.) plus original data.
• Uses augment from broom package
• Idea: “augment” regression with data:

DC.1 %>% augment(windmill) -> DC.1a

• Then use “augmented” dataframe as raw material for plots.

The augmented data

glimpse(DC.1a)

Rows: 25
Columns: 8
$ wind_velocity <dbl> 5.00, 6.00, 3.40, 2.70, 10.00, 9.70,…
$ DC_output     <dbl> 1.582, 1.822, 1.057, 0.500, 2.236, 2…
$ .fitted       <dbl> 1.3366195, 1.5777683, 0.9507813, 0.7…
$ .resid        <dbl> 0.24538052, 0.24423165, 0.10621871, …
$ .hat          <dbl> 0.04834508, 0.04011347, 0.08860703, …
$ .sigma        <dbl> 0.2353240, 0.2354330, 0.2401887, 0.2…
$ .cooksd       <dbl> 0.0288421683, 0.0233028359, 0.010799…
$ .std.resid    <dbl> 1.0655959, 1.0560493, 0.4713466, -1.…

Plot of residuals against fitted values

ggplot(DC.1a, aes(y = .resid, x = .fitted)) + geom_point()

Comments on residual plot

• Residual plot should be a random scatter of points.
• Should be no pattern “left over” after fitting the regression.
• Smooth trend should be more or less straight across at 0.
• Here, have a curved trend on residual plot.
• This means original relationship must have been a curve (as we saw on original scatterplot).
• Possible ways to model a curve:
  • Add a squared term in explanatory variable.
  • Transform response variable (doesn’t work well here).
  • See what science tells you about mathematical form of relationship, and try to apply.

Normal quantile plot of residuals

ggplot(DC.1a, aes(sample = .resid)) + 
  stat_qq() + stat_qq_line()

Comments

• normal quantile plot of residuals is slightly left-skewed
• bigger problem is curve on residuals vs. fitted
• so we should try to model curve.

Parabolas and fitting parabola model

• A parabola has equation \[y = ax^2 + bx + c\] with coefficients \(a, b, c\). About the simplest function that is not a straight line.
• Fit one using lm by adding \(x^2\) to right side of model formula with +:

DC.2 <- lm(DC_output ~ wind_velocity + I(wind_velocity^2),
  data = windmill
)

• The I() necessary because ^ in model formula otherwise means something different (to do with interactions in ANOVA).
• Call it parabola model.

Parabola model output

summary(DC.2)

Call:
lm(formula = DC_output ~ wind_velocity + I(wind_velocity^2), 
    data = windmill)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.26347 -0.02537  0.01264  0.03908  0.19903 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -1.155898   0.174650  -6.618 1.18e-06 ***
wind_velocity       0.722936   0.061425  11.769 5.77e-11 ***
I(wind_velocity^2) -0.038121   0.004797  -7.947 6.59e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1227 on 22 degrees of freedom
Multiple R-squared:  0.9676,    Adjusted R-squared:  0.9646 
F-statistic: 328.3 on 2 and 22 DF,  p-value: < 2.2e-16

Comments on output

• R-squared has gone up a lot, from 0.874 (line) to 0.968 (parabola).
• Coefficient of squared term strongly significant (P-value \(6.59 \times 10^{−8}\)).
• Adding squared term has definitely improved fit of model.
• Parabola model better than linear one.
• But…need to check residuals again.

Augment

DC.2 %>% augment(windmill) -> DC.2a

Residual plot from parabola model

ggplot(DC.2a, aes(y = .resid, x =  .fitted)) +
  geom_point()

Normal quantile plot of residuals

ggplot(DC.2a, aes(sample = .resid)) + stat_qq() + 
  stat_qq_line()

This distribution has long tails, which should worry us at least some.

Make scatterplot with fitted line and curve

ggplot(windmill, aes(y = DC_output, x = wind_velocity)) +
  geom_point() + 
  geom_line(data = DC.1a, aes(y = .fitted), 
            colour = "red") +
  geom_line(data = DC.2a, aes(y = .fitted), 
            colour = "blue") -> g

Comments

• This plots:
  • scatterplot (geom_point);
  • fitted values from straight line (via DC.1a);
  • fitted values from parabola (via DC.2a)
  • these last two joined by lines (with points not shown).
• Trick in the geom_line: use the predictions as the y-values to join by lines (from DC.1a and DC.2a), instead of the original data points. Without the data and aes in the geom_lines, original data points (not predictions) would be joined by lines.

Scatterplot with fitted line and curve

Curve clearly fits better than line.

Another approach to a curve

• There is a problem with parabolas, which we’ll see later.

• Ask engineer, “what should happen as wind velocity increases?”:

  • Upper limit on electricity generated, but otherwise, the larger the wind velocity, the more electricity generated.

• Mathematically, asymptote. Straight lines and parabolas don’t have them, but eg. \(y = 1/x\) does: as \(x\) gets bigger, \(y\) approaches zero without reaching it.

• What happens to \(y = a + b(1/x)\) as \(x\) gets large?

  • \(y\) gets closer and closer to \(a\): that is, \(a\) is asymptote.

• Fit this, call it asymptote model.

• Fitting the model here because we have math to justify it.

  • Alternative, \(y = a + be^{−x}\) , approaches asymptote faster.

Plotting DC output against reciprocal of wind velocity

ggplot(windmill, aes(y = DC_output, 
                     x = 1 / wind_velocity)) +
  geom_point() + geom_smooth(se = F)

This is very close to straight.

Comments

• Define new explanatory variable to be \(1/x\), and predict \(y\) from it.

• \(x\) is velocity, distance over time.

• So \(1/x\) is time over distance. In walking world, if you walk 5 km/h, take 12 minutes to walk 1 km, called your pace.

• in a model formula, / has special meaning (as ^ does)

• so in lm, wrap 1 / wind_velocity in I() (like I(wind_velocity^2)).

Then, run regression like this:

DC.3 <- lm(DC_output ~ I(1/wind_velocity), data = windmill)
summary(DC.3)

Call:
lm(formula = DC_output ~ I(1/wind_velocity), data = windmill)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.20547 -0.04940  0.01100  0.08352  0.12204 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)          2.9789     0.0449   66.34   <2e-16 ***
I(1/wind_velocity)  -6.9345     0.2064  -33.59   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.09417 on 23 degrees of freedom
Multiple R-squared:   0.98, Adjusted R-squared:  0.9792 
F-statistic:  1128 on 1 and 23 DF,  p-value: < 2.2e-16

Comments

• R-squared, 0.98, even higher than for parabola model, 0.968.
• Simpler model, only one explanatory variable vs. 2 for parabola model (wind.velocity and its square).
• 1 / wind_velocity (unsurprisingly) strongly significant.
• Looks good, but check residual plot (over).

Residual plot for asymptote model

DC.3 %>% augment(windmill) -> DC.3a
ggplot(DC.3a, aes(y = .resid, x = .fitted)) + geom_point()

Normal quantile plot of residuals

ggplot(DC.3a, aes(sample = .resid)) + 
  stat_qq() + stat_qq_line()

This is skewed (left), but is not bad (and definitely better than the one for the parabola model).

Comments

• Residual plot not bad. But residuals go up to 0.10 and down to −0.20, suggesting possible skewness (not normal). I think it’s not perfect, but OK overall.
• Next: plot scatterplot with all three fitted lines/curves on it (for comparison).

Making the plot

• rather repetitive, but at least easier to see what is going on:

ggplot(windmill, aes(x = wind_velocity, y = DC_output)) + 
  geom_point() +
  geom_line(data = DC.1a, aes(y = .fitted), colour = "red") +
  geom_line(data = DC.2a, aes(y = .fitted), colour = "darkgreen") +
  geom_line(data = DC.3a, aes(y = .fitted), colour = "blue") -> g3

The plot

Comments

• Predictions from curves are very similar.
• Predictions from asymptote model as good, and from simpler model (one \(x\) not two), so prefer those.
• Go back to asymptote model summary.

Asymptote model summary

tidy(DC.3)

Comments

• Intercept in this model about 3.
• Intercept of asymptote model is the asymptote (upper limit of DC.output).
• Not close to asymptote yet.
• Therefore, from this model, wind could get stronger and would generate appreciably more electricity.
• This is extrapolation! Would like more data from times when wind.velocity higher.
• Slope -7. Why negative?
  • As wind velocity increases, its reciprocal goes down, and DC.output goes up. Check.
• Actual slope number hard to interpret.

Checking back in with research questions

• Is there a relationship between wind speed and current generated?
  • Yes.
• If so, what kind of relationship is it?
  • One with an asymptote.
• Can we model the relationship, in such a way that we can do predictions?
  • Yes, see model DC.3 and plot of fitted curve.
• Good. Job done.

Job done, kinda

• Just because the parabola model and asymptote model agree over the range of the data, doesn’t necessarily mean they agree everywhere.
• Extend range of wind.velocity to 1 to 16 (steps of 0.5), and predict DC.output according to the two models:

wv <- seq(1, 16, 0.5)
wv

 [1]  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0
[12]  6.5  7.0  7.5  8.0  8.5  9.0  9.5 10.0 10.5 11.0 11.5
[23] 12.0 12.5 13.0 13.5 14.0 14.5 15.0 15.5 16.0

• R has predict, which requires values to predict for, as data frame.
• This has to contain values, with matching names, for all explanatory variables in regression(s).

Setting up data frame to predict from

• Linear model had just wind_velocity.
• Parabola model had that as well (squared one will be calculated)
• Asymptote model had just reciprocal of velocity, but worked it out from wind_velocity.
• So create data frame called wv_new with values in:

wv_new <- tibble(wind_velocity = wv)

wv_new

wv_new

Doing predictions, one for each model

wv_new %>% 
  mutate(linear = predict(DC.1, wv_new),
         parabola = predict(DC.2, wv_new),
         asymptote = predict(DC.3, wv_new)) -> my_fits

my_fits

my_fits

Two ways to make a plot of these:

• the “Excel method”, plotting the columns one at a time
• the “ggplot method”, to plot one thing labelled different ways

The “Excel method”

ggplot(my_fits, aes(x = wind_velocity, y = linear)) +
  geom_line(colour = "red") +
  geom_line(aes(y = parabola), colour = "darkgreen") +
  geom_line(aes(y = asymptote), colour = "blue") -> g_excel

The plot

Comments

• This coding is easy enough to follow, but rather repetitive
• When you find yourself repeating things, time to wonder whether there is a better way
• ggplot likes having the x and y for a plot in single columns, with additional columns labelling things (like the model, here)
• This is what pivot_longer can help us with.

The “ggplot method”

my_fits %>%
    pivot_longer(linear:asymptote,
      names_to="model", values_to="fit") -> my_fits_longer
my_fits_longer

Now the plotting is much simpler

ggplot(my_fits_longer, aes(y = fit, 
                           x = wind_velocity,
                           colour = model)) +
  geom_line() -> ggg

The plot so far

Showing where we had data

• Would be nice to see where our original data were.
• The observed wind velocities were in this range:

(vels <- range(windmill$wind_velocity))

[1]  2.45 10.20

• DC.output between 0 and 3 from asymptote model. Add rectangle to graph around where the data were:

ggg + geom_rect(
  xmin = vels[1], xmax = vels[2], ymin = 0, ymax = 3,
  alpha = 0, colour = "black") -> gr

The plot now

Comments (1)

• Over range of data, two models agree with each other well.
• Outside range of data, they disagree violently!
• For larger wind.velocity, asymptote model behaves reasonably, parabola model does not.
• What happens as wind.velocity goes to zero? Should find DC.output goes to zero as well. Does it?

Comments (2)

• For parabola model:

tidy(DC.2)

• Nope, goes to −1.16 (intercept), actually significantly different from zero.

Comments (3): asymptote model

tidy(DC.3)

• As wind.velocity heads to 0, its reciprocal heads to \(+\infty\), so DC.output heads to \(−\infty\)!
• Also need more data for small wind.velocity to understand relationship. (Is there a lower asymptote?)
• Best we can do now is to predict DC.output to be zero for small wind.velocity.
• Assumes a “threshold” wind velocity below which no electricity generated at all.

Summary

• Often, in data analysis, there is no completely satisfactory conclusion, as here.
• Have to settle for model that works OK, with restrictions.
• Always something else you can try.
• At some point you have to say “I stop.”