Assignment 5 - 2011 - Solution

Due date(s):	16 February 2011
(PDF)	Assignment questions
(PDF)	Assignment solutions

<rst> <rst-options: 'toc' = False/> <rst-options: 'reset-figures' = True/>

Assignment objectives

=========

.. rubric:: Assignment objective: explore and understand least squares models a bit more.

.. rubric:: Brief solutions (i.e. using bullet points and short explanations) are provided for questions 1 and 2, but they are complete. You may answer questions in the midterm and final exam at this level.

Question 1 [2]

=====

Use the `bioreactor data <http://openmv.net/info/bioreactor-yields>`_, which shows the percentage yield from the reactor when running various experiments where temperature was varied, impeller speed and the presence/absence of baffles were adjusted.

1. . Build a linear model that uses the reactor temperature to predict the yield. Interpret the slope and intercept term.

1. . Build a linear model that uses the impeller speed to predict yield. Interpret the slope and intercept term.

1. . Build a linear model that uses the presence (represent it as 1) or absence (represent it as 0) of baffles to predict yield. Interpret the slope and intercept term.

*Note*: if you use R it will automatically convert the ``baffles`` variable to 1's and 0's for you. If you wanted to make the conversion yourself, to verify what R does behind the scenes, try this:

.. code-block:: s

# Read in the data frame bio <- read.csv('http://datasets.connectmv.com/file/bioreactor-yields.csv')

# Force the baffles variables to 0's and 1's bio$baffles <- as.numeric(bio$baffles) - 1

1. . Which variable(s) would you change to boost the batch yield, at the lowest cost of implementation?

1. . Use the ``plot(bio)`` function in R, where ``bio`` is the data frame you loaded using the ``read.csv(...)`` function. R notices that ``bio`` is not a single variable, but a group of variables, i.e. a data frame, so it plots what is called a *scatterplot matrix* instead. Describe how the scatterplot matrix agrees with your interpretation of the slopes in parts 1, 2 and 3 of this question.

Solution

---

The R code (below) was used to answer all questions.

1. .

* The model is: :math:`\hat{y} = 102.5 - 0.69T`, where :math:`T` is tank temperature. * Intercept = :math:`102.5` % points is the yield when operating at 0 :math:`^\circ \text{C}`. Obviously not a useful interpretation, because data have not been collected in a range that spans, or is even close to 0 :math:`^\circ \text{C}`. It is likely that this bioreactor system won't yield any product under such cold conditions. Further, a yield greater than 100% is not realizable. * Slope = -0.69 :math:`\frac{[\%]}{[^\circ \text{C}]}`, indicating the yield decreases, on average, by about 0.7 units for every degree increase in tank temperature.

1. .

* The model is: :math:`\hat{y} = -20.3 + 0.016S`, where :math:`S` is impeller speed. * Intercept = :math:`-20.3` % points is the yield when operating no agitation. Again, obviously not a useful interpretation, because the data have not been collected under these conditions, and yield can't be a negative quantity. * Slope = 0.016 :math:`\frac{[\%]}{[\text{RPM}]}`, indicating the yield increases, on average, by about 1.6 percentage points per 100 RPM increase.

1. .

* The model is: :math:`\hat{y} = 54.9 - 16.7B`, where :math:`B` is 1 if baffles are present and :math:`B=0` with no baffles. * Intercept = :math:`54.9` % points yield is the yield when operating with no baffles (it is in fact the average yield of all the rows that have "No" as their baffle value). * Slope = -16.7 %, indicating the presence of baffles decreases the yield, on average, by about 16.7 percentage points.

1. . This is of course an open-ended, and case specific. Some factors you would include are:

* Remove the baffles, but take into account the cost of doing so. Perhaps it takes a long time (expense) to remove them, especially if the reactor is used to produce other products that do require the baffles.

* Operate at lower temperatures. The energy costs of cooling the reactor would factor into this.

* Operate at higher speeds and take that cost into account. Notice however there is one observation at 4900 RPM that seems unusual: was that due to the presence of baffles, or due to temperature in that run? We'll look into this issue with multiple linear regression later on.

.. note::

Please note that our calculations above are not the true effect of each of the variables (temperature, speed and baffles) on yield. Our calculations assume that there is no interaction between temperature, speed and baffles, and that each effect operates independent of the others. That's not necessarily true. After the midterm break we will learn how to "control for the effects" of other variables and build a single model with all variables.

1. . The scatterplot matrix, shown below, agrees with our interpretation. This is an information rich visualization that gives us a feel for the multivariate relationships and really summarizes all the variables well (especially the last row of plots).

* The yield-temperature relationship is negative, as expected. * The yield-speed relationship is positive, as expected. * The yield-baffles relationship is negative, as expected. * We can't tell anything about the yield-duration relationship, as it doesn't vary in the data we have (there could/should be a relationship, but we can't tell).

.. literalinclude:: ../che4c3/Assignments/Assignment-5/code/bioreactor-regression-assignment.R :language: s

.. figure:: ../che4c3/Assignments/Assignment-5/images/bioreactor-scatterplot-matrix-assignment.png :alt: images/math :width: 550px :align: center

Question 2 [2]

=====

Use the `gas furnace data <http://datasets.connectmv.com/info/gas-furnace>`_ from the website to answer these questions. The data represent the gas flow rate (centered) from a process and the corresponding CO\ :sub:`2` measurement.

1. . Make a scatter plot of the data to visualize the relationship between the variables. How would you characterize the relationship?

1. . Calculate the variance for both variables, the covariance between the two variables, and the correlation between them, :math:`r(x,y)`. Interpret the correlation value; i.e. do you consider this a strong correlation?

1. . Now calculate a least squares model relating the gas flow rate as the :math:`x` variable to the CO\ :sub:`2` measurement as the :math:`y`-variable. Report the intercept and slope from this model.

1. . Report the :math:`R^2` from the regression model. Compare the squared value of :math:`r(x,y)` to :math:`R^2`. What do you notice? Now reinterpret what the correlation value means (i.e. compare this interpretation to your answer in part 2).

1. . **600-level**: Switch :math:`x` and :math:`y` around and rebuild your least squares model. Compare the new :math:`R^2` to the previous model's :math:`R^2`. Is this result surprising? How do interpret this?

Solution

---

1. . Relationship: the data are negatively correlated.

.. figure:: ../che4c3/Assignments/Assignment-5/images/CO2-gas-rate-raw-data-assignment5.png :alt: ../che4c3/Assignments/Assignment-5/code/gas-furnace-question-assignment2.R :width: 450px :align: center

I've chosen to use the ``sp`` or ``scatterplot`` function from the ``car`` library. It shows the scatterplot smoother (a.k.a. loess line) as solid red, the spread around the smoother (dashed red), the least squares regression line (black) and boxplots for each axis.

This is a great example of an information-rich visualization: packing the maximum amount of information into a small space. This plot answers so many questions we might have about the data.

1. . The ``cov(...)`` command supplies the variance and covariance, and the ``cor(...)`` command gives the correlation.

* Variance of input gas flow rate = 1.15 [gas flow units] :math:`^2` * Variance of CO\ :sub:`2` = 10.3 [CO\ :sub:`2` units] :math:`^2` * Covariance between input gas flow and CO\ :sub:`2` = -1.66 [gas flow units][CO\ :sub:`2` units] * Correlation = -0.48, i.e. around -0.5.

From my experience with data, I personally would interpret this as a reasonably strong correlation. There is reasonably strong linear behaviour in the data cloud shown above, enough of a relationship to confidently say that "the CO\ :sub:`2` output does decrease at higher gas flow rates".

1. . From the R model output:

	*	intercept is -1.44 units of CO\ :sub:`2` 

* slope is 53.4 :math:`\frac{[\text{units of CO}_2]}{[\text{units of gas flow}]}`

1. .

* From the R model output: :math:`R^2 = 0.2347` * From earlier, the squared correlation is :math:`(-0.484)^2 = 0.2347`, the same value. * Correlation can be interpreted as the square root of the :math:`R^2` value when regressing :math:`y` on :math:`x` (i.e. fitting a linear model to :math:`y` using :math:`x` as the input). * Most novices would be misled and consider an :math:`R^2` value of 0.23 quite low. But notice that there is a repeatable and consistent negative linear relationship between :math:`x` and :math:`y` in this data.

1. . This shows the interesting result that when regressing :math:`x` on :math:`y` (instead of the usual regression of :math:`y` on :math:`x`), that we get the same :math:`R^2` value. Note however that the *intercept* and *slope* are different between the two regressions.

This also calls into question the interpretation of the :math:`R^2` value in regression. :math:`R^2` is just the square of the correlation coefficient. Recall from class the slide on the `Wikipedia examples of correlation <http://en.wikipedia.org/wiki/File:Correlation_examples.png>`_: there were examples where :math:`r(x,y) = \sqrt{R^2}` was zero, but still a strong *relationship* existing in the data. So we should interpret :math:`R^2` as a measure only of the *linear relationship* between two variables. And bear its quadratic nature in mind - interpreting the correlation is actually easier, and more "linear", in that a 0.2 improvement in correlation means the same thing when going from :math:`r=0.2` to 0.4, as it does when going from :math:`r=0.7` to 0.9 (not so for :math:`R^2`).

.. literalinclude:: ../che4c3/Assignments/Assignment-5/code/gas-furnace-question-assignment5.R :language: s

Question 3 [1.5]

====

A new type of `thermocouple <http://en.wikipedia.org/wiki/Thermocouple>`_ is being investigated by your group. These devices produce an *almost* linear voltage (millivolt) response at different temperatures. In practice though it is used the other way around: use the millivolt reading to predict the temperature. The process of fitting this linear model is called *calibration*.

1. . Use the following data to calibrate a linear model:

================= ==== ==== ==== ==== ==== ==== ==== ==== ==== ==== Temperature [K] 273 293 313 333 353 373 393 413 433 453 ----------------- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- Reading [mV] 0.01 0.12 0.24 0.38 0.51 0.67 0.84 1.01 1.15 1.31 ================= ==== ==== ==== ==== ==== ==== ==== ==== ==== ====

Show the linear model and provide the predicted temperature when reading 1.00 mV.

1. . Are you satisfied with this model, based on the coefficient of determination (:math:`R^2`) value?

1. . What is the model's standard error? Now, are you satisfied with the model's prediction ability, given that temperatures can usually be recorded to an accuracy of :math:`\pm 0.5` K with most inexpensive thermocouples.

1. . What is your (revised) conclusion now about the usefulness of the :math:`R^2` value?

• • Note**: This example explains why I don't use the terminology of *independent* and *dependent* variables in this course. Here the temperature truly is the independent variable, because it causes the voltage difference that we measure. But the voltage reading is the independent variable in the least squares model. The word *independent* is being used in two different senses (its English meaning *vs* its mathematical meaning), and this can be misleading.

Solution

---

1. . The linear model is used to predict temperature given the reading in millivolts. The reason is that in modelling, in general, we specify as :math:`x` the variable(s) we always have available, while :math:`y` is the variable we would like to predict from the :math:`x`.

The model has the form: :math:`T = b_0 + b_1V`, where :math:`T` is temperature and :math:`V` is the voltage reading. Coefficients in the linear model are:

.. math::

T = 278.6 + 135.3 V

implies that recording an increase in 0.1 mV means, on average, the temperature has increased by 13.5 K in the system.

The temperature prediction at 1.00 mV would be 413.9 K.

.. figure:: ../che4c3/Assignments/Assignment-5/images/voltage-linear-model.png :width: 450px :align: center

The following code was used to fit the model and draw the plot.

.. twocolumncode:: :code1: ../che4c3/Assignments/Assignment-5/code/voltage_linear_model.m :language1: matlab :header1: MATLAB code :code2: ../che4c3/Assignments/Assignment-5/code/voltage_linear_model.py :language2: python :header2: Python code

If you used ``R`` to fit the model, you would written something like this::

> V <- c(0.01, 0.12, 0.24, 0.38, 0.51, 0.67, 0.84, 1.01, 1.15, 1.31) > T <- c(273, 293, 313, 333, 353, 373, 393, 413, 433, 453) > model <- lm(T ~ V) > summary(model)

Call: lm(formula = T ~ V)

Residuals: Min 1Q Median 3Q Max -6.9272 -2.1212 -0.1954 2.7480 5.4239

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 278.574 2.204 126.39 1.72e-14 *** V 135.298 2.922 46.30 5.23e-11 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.916 on 8 degrees of freedom Multiple R-squared: 0.9963, Adjusted R-squared: 0.9958 F-statistic: 2144 on 1 and 8 DF, p-value: 5.229e-11

1. . The :math:`R^2` value from this linear fit is :math:`R^2 = 0.996`, which being so close to 1.0, implies the linear relationship in the data is strong (the linear model fits the data very well) - that's all.

One cannot be satisfied with only an :math:`R^2` value: it has nothing to do with whether the model's prediction accuracy is any good. So we can't tell anything from this number.

1. . The model's standard error is 3.9 K. If we assume the prediction error is normally distributed around the linear fit, this corresponds to one standard deviation. So 95% of our prediction error lies roughly within a range of :math:`\pm 2\times 3.92` or :math:`\pm 7.8` K. These are the dashed red lines drawn on the figure. (Please note: the true error intervals are not parallel to the regression line, they are curved; however the :math:`\pm 2S_E` limits are a good-enough approximation for most engineering applications.

This prediction ability of :math:`\pm 8` K is probably not satisfying for most engineering applications, since we can predict temperatures far more accurately, over the range from 273K to 453K, using off-the-shelf commercial thermocouples.

1. . The purpose of this question is to mainly point out the misleading nature of :math:`R^2` - this value looks really good: 99.6%, yet the actual purpose of the model, the ability to predict temperature from the millivolt reading, has no relationship at all to this :math:`R^2` value.

.. sd(T) = 60.5 .. diff(range(T)) = 180 .. baseline_ratio = 60/180 = 0.3333 .. SE = 3.9 .. model_ratio = 3.9/180 = 0.0216 .. ratio = 1 - 0.0216/0.3333 = 0.935: seems pretty good

Question 4 [1]

=====

1. . Use the linear model you derived in Question 2, where you used the gas flow rate to predict the CO\ :sub:`2` measurement, and construct the analysis of variance table (ANOVA) for the dataset. Use your ANOVA table to reproduce the residual standard error, :math:`S_E` value, that you get from the R software output.

Go through the `R tutorial <http://connectmv.com/tutorials/r-tutorial/>`_ to learn how to efficiently obtain the residuals and predicted values from a linear model object.

1. . Also for linear model in Question 2, verify whether the residuals are normally distributed.

1. . Use the linear model you derived in Question 3, where you used the voltage measurement to predict the temperature, and construct the analysis of variance table (ANOVA) for that dataset. Use your ANOVA table to reproduce the residual standard error, :math:`S_E` value, that you get from the R software output.

.. raw:: latex

\newpage

Solution

---

1. . The ANOVA table values were calculated in the code solutions for question 2:

=================== ========================================= =================== ============== ======================================== Type of variance Distance Degrees of freedom SSQ Mean square =================== ========================================= =================== ============== ======================================== Regression :math:`\hat{y}_i - \overline{\mathrm{y}}` :math:`k-2` 709.9 354.9 ------------------- ----------------------------------------- ------------------- -------------- ---------------------------------------- Error :math:`y_i - \hat{y}_i` :math:`n-k` 2314.9 7.87 ------------------- ----------------------------------------- ------------------- -------------- ---------------------------------------- Total :math:`y_i - \overline{\mathrm{y}}` :math:`n` 3024.8 10.2 =================== ========================================= =================== ============== ========================================

The residual standard error, or just standard error, :math:`S_E = \sqrt{\frac{2314.9}{296-2}} = 2.8` %CO\ :sub:`2`, which agrees with the value from R.

1. . These residuals are normally distributed, as verified in the following qq plot:

.. figure:: ../che4c3/Assignments/Assignment-5/images/residuals-gas-furnace-question-assignment.png :alt: ../che4c3/Assignments/Assignment-5/code/gas-furnace-question-assignment2.R :width: 450px :align: center

As mentioned in the ``help(qqPlot)`` output, the dashed red line is the confidence envelope at the 95% level. The single point just outside the confidence envelope is not going to have any practical effect on our assumption of normality. We expect 1 point in 20 to lie outside the limits.

We will discuss the meaning of "studentized residuals", on the :math:`y`-axis, after the midterm.

1. . For the question 3 data set:

=================== ========================================= =================== ============== ======================================== Type of variance Distance Degrees of freedom SSQ Mean square =================== ========================================= =================== ============== ======================================== Regression :math:`\hat{y}_i - \overline{\mathrm{y}}` :math:`k-2` 32877 16438 ------------------- ----------------------------------------- ------------------- -------------- ---------------------------------------- Error :math:`y_i - \hat{y}_i` :math:`n-k` 122.7 15.3 ------------------- ----------------------------------------- ------------------- -------------- ---------------------------------------- Total :math:`y_i - \overline{\mathrm{y}}` :math:`n` 33000 3300 =================== ========================================= =================== ============== ========================================

The residual standard error, or just standard error, :math:`S_E = \sqrt{\frac{122.7}{10-2}} = 3.9` K, which agrees with the value from R.

.. raw:: latex

\vspace{0.5cm} \hrule \begin{center}END\end{center}

</rst>