class: split-70 with-border hide-slide-number bg-brand-red background-image: url("images/USydLogo-black.svg") background-size: 200px background-position: 2% 90% .column.white[.content[ <br><br><br> # Polynomial Regression ## .black[STAT3022 Applied Linear Models Lecture 12] <br><br><br> ### .black[2020/02/17] ]] .column.bg-brand-charcoal[.content.white[ ## Today 1. Polynomial regression 1. Orthogonal polynomials 1. The poly, vcov and cov2cor command ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Polynomial functions: mix & match ]] .row[.split-two[ .row[.split-five[ .column[.content.font_large[ A) `$$y = x^2$$` ]] .column.bg-brand-gray[.content.font_large[ B) `$$y = 3 + x^3$$` ]] .column[.content[ .font_large[ C)] .font_medium[ `$$y = x^6 - 9x^3 + 8$$` ]]] .column.bg-brand-gray[.content.font_large[ D) `$$y = 3x + 2$$` ]] .column[.content.font_large[ E) `$$y = x^6 - 8x^3$$` ]] ]] .row[.split-five[ .column.bg-brand-gray[.content.font_large[ I) <img src="lecture12_2020JC_files/figure-html/unnamed-chunk-3-1.svg" style="display: block; margin: auto;" /> ]] .column[.content.font_large[ II) <img src="lecture12_2020JC_files/figure-html/unnamed-chunk-4-1.svg" style="display: block; margin: auto;" /> ]] .column.bg-brand-gray[.content.font_large[ III) <img src="lecture12_2020JC_files/figure-html/unnamed-chunk-5-1.svg" style="display: block; margin: auto;" /> ]] .column[.content.font_large[ IV) <img src="lecture12_2020JC_files/figure-html/unnamed-chunk-6-1.svg" style="display: block; margin: auto;" /> ]] .column.bg-brand-gray[.content.font_large[ V) <img src="lecture12_2020JC_files/figure-html/unnamed-chunk-7-1.svg" style="display: block; margin: auto;" /> ]] ]]]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # The polynomial regression model ]] .row[.split-60[ .column[.content[ * A polynomial linear regression model of degree `\(k\in\Bbb{N}\)` has the form `$$Y_i = \beta_0 + \beta_1 x_{i} + \beta_2 x_i^2 + \ldots + \beta_k x_{i}^k + \epsilon_i, ~~~~ (i=1,\ldots,n),$$` where `\(\epsilon_1, \ldots \epsilon_n\)` are random errors satisfying the usual assumptions (A1) - (A4). * Typically statisticians *aim for low order polynomial models* - having degree `\(k\)` larger than 5 or 6 is unusual. ### Include lower order terms! .blockquote[ If you use a polynomial model of degree `\(k\)`, then you should include all terms of lower order (unless the circumstances are exceptional). ] ]] .column.bg-brand-gray[.content[ ### Example: Cubic polynomials * A cubic polynomial should include linear and quadratic terms: `$$E[Y_i | x_{i}] = \beta_0 + \beta_1 x_{i} + \beta_2 x_i^2 + \beta_3 x_{i}^3$$` * We should *not* omit, say, the quadratic term: `$$E[Y_i | x_{i}] = \beta_0 + \beta_1 x_{i} + \beta_3 x_{i}^3$$` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Child lung function data ]] .row[.split-30[ .column[.content[ * FEV (forced expiratory volume) is a measure of lung function. * The data include determinations of FEV on 318 female children who were seen in a childhood respiratory disease study in Massachusetts in the U.S. * Child's age (in years) also recorded. * Of interest is to model FEV by age. .bottom_abs.width100.font_small[ Reference: Tager et al. (1979). *American J of Epidemiology*, **110**, 15-26. ]]] .column.bg-brand-gray[.content[ ```r dat <- read.table(file = "data/fev.txt", header = T) skimr::skim(dat) ``` ``` Skim summary statistics n obs: 318 n variables: 2 -- Variable type:integer --------------------------------------------------------------------------------------------------- variable missing n mean sd p0 p25 p50 p75 p100 hist Age 0 318 9.84 2.93 3 8 10 12 19 <U+2582><U+2583><U+2587><U+2587><U+2585><U+2582><U+2581><U+2581> -- Variable type:numeric --------------------------------------------------------------------------------------------------- variable missing n mean sd p0 p25 p50 p75 p100 hist FEV 0 318 2.45 0.65 0.79 1.95 2.49 2.99 3.83 <U+2581><U+2583><U+2585><U+2585><U+2587><U+2586><U+2585><U+2582> ``` ```r cor(dat$FEV, dat$Age) ``` ``` [1] 0.7390163 ``` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Lung data: Scatterplot ]] .row[.split-two[ .column.bg-brand-gray[.content[ <br><br> <img src="lecture12_2020JC_files/figure-html/unnamed-chunk-9-1.svg" style="display: block; margin: auto;" /> <br> Does the linear model look like a good fit to the data? ]] .column[.content[ <br> <img src="lecture12_2020JC_files/figure-html/unnamed-chunk-10-1.png" style="display: block; margin: auto;" /> ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Lung data: Degree 6 polynomial regression model ]] .row[.split-60[ .column[.content[ ```r dat <- dat %>% mutate(Age2=Age^2, Age3=Age^3, Age4=Age^4, Age5=Age^5, Age6=Age^6) MP6 <- lm(FEV ~ Age + Age2 + Age3 + Age4 + Age5 + Age6, data=dat) broom::tidy(MP6) ``` ``` term estimate std.error statistic p.value 1 (Intercept) -1.519931e+00 5.052333e+00 -0.3008376 0.7637394 2 Age 2.078781e+00 3.532837e+00 0.5884171 0.5566795 3 Age2 -6.468470e-01 9.734634e-01 -0.6644800 0.5068755 4 Age3 1.033263e-01 1.358506e-01 0.7605877 0.4474796 5 Age4 -8.144987e-03 1.016624e-02 -0.8011797 0.4236393 6 Age5 3.056270e-04 3.882215e-04 0.7872490 0.4317356 7 Age6 -4.359556e-06 5.930743e-06 -0.7350776 0.4628462 ``` ]] .column.bg-brand-gray[.content[ * We started by trying a polynomial model of degree 6, `MP6`. * The `\(x^6\)` term can be removed (coefficient of `Age6` not significant). * Note that the linear term `Age` has largest p-value but focus centers on the highest order term first. ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Lung data: Simplify further ]] .row[.split-two[ .column[.content[ ```r MP5 <- update(MP6, . ~ . - Age6) print(broom::tidy(MP5), digits = 4) ``` ``` term estimate std.error statistic p.value 1 (Intercept) 1.750e+00 2.393e+00 0.7312 0.4652 2 Age -3.209e-01 1.349e+00 -0.2378 0.8122 3 Age2 3.714e-02 2.858e-01 0.1299 0.8967 4 Age3 5.707e-03 2.860e-02 0.1995 0.8420 5 Age4 -7.392e-04 1.359e-03 -0.5438 0.5870 6 Age5 2.083e-05 2.467e-05 0.8444 0.3991 ``` * The 5th order term can be removed (coefficient of `Age5` not significant). * We try a polynomial model of degree 4, `MP4`: ]] .column.bg-brand-gray[.content[ ```r MP4 <- update(MP5, . ~ . - Age5) print(broom::tidy(MP4), digits=2) ``` ``` term estimate std.error statistic p.value 1 (Intercept) 3.56859 1.0434 3.4 7.1e-04 2 Age -1.39412 0.4529 -3.1 2.3e-03 3 Age2 0.27128 0.0692 3.9 1.1e-04 4 Age3 -0.01815 0.0044 -4.1 5.4e-05 5 Age4 0.00041 0.0001 4.0 7.1e-05 ``` ```r outMP4 <- broom::tidy(MP4) ``` * The `\(x^4\)` term is statistically significant, so this is our preferred model and the fitted model equation is: `$$\hat{E}[\mbox{FEV}] = 3.57 - 1.39~\mbox{Age} + 0.27~\mbox{Age}^2 - 0.018~\mbox{Age}^3$$` `$$\qquad + 0.00041~\mbox{Age}^4.$$` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Lung data: Remarks on the final model ]] .row[.split-two[ .column[.content[ ```r data.frame(x=c(min(dat$Age), max(dat$Age))) %>% ggplot(aes(x)) + stat_function(fun=function(x) outMP4$estimate[1] + outMP4$estimate[2] * x + outMP4$estimate[3] * x ^ 2 + outMP4$estimate[4] * x ^ 3 + outMP4$estimate[5] * x ^ 4 ) + geom_point(data=dat, aes(Age, FEV)) + theme_bw(base_size=18) ``` <img src="lecture12_2020JC_files/figure-html/unnamed-chunk-15-1.png" style="display: block; margin: auto;" /> ]] .column.bg-brand-gray[.content[ * Looking at the plot, the fitted model looks adequate for ages in the range 5 to 15 years. * At the extreme ends the fitted curve exhibits spurious upward curvature but there is insufficient data. * Therefore, the form of the fitted polynomial is an artefact rather than a characteristic of the data. * Models are unstable with coefficients change greatly after dropping each higher order before a suitable functional form is found. Eg Coefficient for .brand-blue[Age] in MP6 = 2.08 and in MP5 = -0.321. ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Collinearity ]] .row[.split-30[ .column[.content[ * There is typically noticeable .brand-blue[collinearity] between the powers of any explanatory variable. * Collinearity in polynomial regression can lead to .brand-blue[numerical instability] in statistical software packages. * Problems of collinearity can be addressed by using .brand-blue[orthogonal polynomials]. ]] .column.bg-brand-gray[.content[ ### Variance-covariance matrix of regression parameter esimates For the linear model `$$\boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\epsilon},\quad \boldsymbol{\epsilon}\sim N(\boldsymbol{0}, \sigma^2\boldsymbol{I}_n).$$` Recall `\(\widehat{\text{V}}\text{ar}(\hat{\boldsymbol{\beta}})= \hat{\sigma}^2(\mathbf{X}^\top\mathbf{X})^{-1}\)` and this can be calculated in R by using `vcov`. ```r round(Sigma <- vcov(MP6), 3) ``` ``` (Intercept) Age Age2 Age3 Age4 Age5 Age6 (Intercept) 25.526 -17.732 4.807 -0.656 0.048 -0.002 0 Age -17.732 12.481 -3.422 0.471 -0.035 0.001 0 Age2 4.807 -3.422 0.948 -0.132 0.010 0.000 0 Age3 -0.656 0.471 -0.132 0.018 -0.001 0.000 0 Age4 0.048 -0.035 0.010 -0.001 0.000 0.000 0 Age5 -0.002 0.001 0.000 0.000 0.000 0.000 0 Age6 0.000 0.000 0.000 0.000 0.000 0.000 0 ``` ```r D <- diag(1 / sqrt(diag(Sigma))) #diag matrix with diag element 1/se. ``` ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ ## Correlation matrix for regression parameter estimates ]] .row[.split-60[ .column[.content[ ```r round(cov2cor(Sigma), 2) ``` ``` (Intercept) Age Age2 Age3 Age4 Age5 Age6 (Intercept) 1.00 -0.99 0.98 -0.96 0.93 -0.91 0.88 Age -0.99 1.00 -0.99 0.98 -0.97 0.95 -0.92 Age2 0.98 -0.99 1.00 -1.00 0.99 -0.97 0.96 Age3 -0.96 0.98 -1.00 1.00 -1.00 0.99 -0.98 Age4 0.93 -0.97 0.99 -1.00 1.00 -1.00 0.99 Age5 -0.91 0.95 -0.97 0.99 -1.00 1.00 -1.00 Age6 0.88 -0.92 0.96 -0.98 0.99 -1.00 1.00 ``` which is the same as below. All have high correlation. ```r round(D %*% Sigma %*% D, 2) ``` ``` [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1.00 -0.99 0.98 -0.96 0.93 -0.91 0.88 [2,] -0.99 1.00 -0.99 0.98 -0.97 0.95 -0.92 [3,] 0.98 -0.99 1.00 -1.00 0.99 -0.97 0.96 [4,] -0.96 0.98 -1.00 1.00 -1.00 0.99 -0.98 [5,] 0.93 -0.97 0.99 -1.00 1.00 -1.00 0.99 [6,] -0.91 0.95 -0.97 0.99 -1.00 1.00 -1.00 [7,] 0.88 -0.92 0.96 -0.98 0.99 -1.00 1.00 ``` ]] .column.bg-brand-gray[.content[ $$ \widehat{\text{C}}\text{orr}(\hat{\boldsymbol{\beta}}) = \text{Diag}(\hat \sigma_i^{-1}) \widehat{\text{C}}\text{ov}(\hat{\boldsymbol{\beta}})\text{Diag}(\hat \sigma_i^{-1})$$ where `\(\text{Diag}(\hat \sigma_i^{-1})\)` is a diagonal matrix with diagonal entities `\(\hat \sigma_i^{-1}\)` and `\(\widehat{\text{C}}\text{orr}(\beta_i,\beta_j) = \widehat{\text{C}}\text{ov}(\beta_i,\beta_j)/\hat \sigma_i \hat \sigma_j\)`. ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Orthogonal polynomials ]] .row[.split-50[ .column[.content[ * Let `\(P_r(x)\)` denote a polynomial of degree `\(r\)` in `\(x\)` and `\(x\)` takes sample values `\(x_1, x_2, \ldots, x_n\)`. * Two polynomials `\(P_r\)` and `\(P_s\)` are *orthogonal* if `$$\sum_{i=1}^n P_r(x_i) P_s(x_i) = 0~~~~~~(r \ne s)$$` * We can fit a polynomial of order `\(k\)` to the data using orthogonal polynomials `\(P_r(x), r=1,\dots,k\)`: `$$E[Y_i] = \alpha_0 + \alpha_1 P_1(x_i) + \alpha_2 P_2(x_i) + \ldots + \alpha_k P_k(x_i).$$` This will be completely equivalent to fitting `$$E[Y_i] = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \ldots + \beta_k x_i^k$$` with the `\(\beta\)`'s being a simple transformation of the `\(\alpha\)`'s. ]] .column.bg-brand-gray[.content[ ### Key property of orthogonal polynomials * Fitting a polynomial regression using orthogonal polynomials is very convenient because the resulting .brand-blue[ estimated regression coefficients are uncorrelated]. * As a result, we can interpret `\(t\)`-test statistics and `\(p\)`-values in the summary table independently, ie .brand-blue[p-values will not change much when we remove other terms from the model]. Good for variable selection! ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Example: Simple linear regression ]] .row[.split-50[ .column[.content[ * Let `\(P_0=1\)` and `\(P_1(x_i) = x_i - \bar x\)`. * The constant polynomial `\(P_0\)` is orthogonal to `\(P_1\)`, since `\(\sum_{i=1}^n P_0(x_i) P_1(x_i) = \sum_{i=1}^n (x_i - \bar x) = 0\)`. * We can fit a simple linear regression model as `$$E[Y_i] = \alpha_0 + \alpha_1 P_1(x_i) = \alpha_0 + \alpha_1 (x_i - \bar x)$$` * This is equivalent to fitting the usual model `$$E[Y_i] = \beta_0 + \beta_1 x_i,$$` where `\(\beta_1 = \alpha_1\)` and `\(\beta_0 = \alpha_0 - \alpha_1 \bar x\)`. * The term `\(\alpha_0\)` represents the mean value of response when `\(x\)` is at its sample mean `\(\bar x\)`. * Estimates `\(\hat \alpha_0\)` and `\(\hat \alpha_1\)` is .brand-blue[uncorrelated] because the .brand-blue[centred] variables `\(z_i = x_i -\bar x\)` have `\(\bar z = 0.\)` ]] .column.bg-brand-gray[.content[ ### Orthogonal polynomials in R Orthogonal polynomials can be fitted in R with the `poly(x,degree=k)` command. ```r MOP6 <- lm(FEV ~ poly(Age, degree = 6), data = dat) print(broom::tidy(MOP6), digits=2) ``` ``` term estimate std.error statistic p.value 1 (Intercept) 2.45 0.022 112.03 2.0e-253 2 poly(Age, degree = 6)1 8.50 0.390 21.78 1.7e-64 3 poly(Age, degree = 6)2 -3.14 0.390 -8.05 1.8e-14 4 poly(Age, degree = 6)3 -0.36 0.390 -0.91 3.6e-01 5 poly(Age, degree = 6)4 1.57 0.390 4.02 7.2e-05 6 poly(Age, degree = 6)5 0.33 0.390 0.84 4.0e-01 7 poly(Age, degree = 6)6 -0.29 0.390 -0.74 4.6e-01 ``` ```r round(broom::tidy(MOP6)$p.value, 3) ``` ``` [1] 0.000 0.000 0.000 0.363 0.000 0.399 0.463 ``` Hence the order of orthogonal polynomial is also 4. ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Lung data: Correlation matrix ]] .row[.split-50[ .column[.content[ ```r outMOP6 <- broom::tidy(MOP6) outMOP6[, c("estimate", "std.error")] ``` ``` estimate std.error 1 2.4511698 0.02188049 2 8.4964720 0.39018494 3 -3.1402745 0.39018494 4 -0.3552419 0.39018494 5 1.5699927 0.39018494 6 0.3292101 0.39018494 7 -0.2868162 0.39018494 ``` * Regression coefficients for `MOP6` are uncorrelated! * Transformed columns of `\(\boldsymbol{X}\)` are uncorrelated and of .brand-blue[equal scale]. Standard errors of `\(\beta_j\)` are equal. * Model fit using `\(R^2\)` for `MP6` and `MOP6` are the same . ```r c(summary(MP6)$r.squared,summary(MOP6)$r.squared) ``` ``` [1] 0.6417945 0.6417945 ``` ]] .column.bg-brand-gray[.content[ ```r round(unname(vcov(MOP6)), 2) #cov matrix of beta ``` ``` [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0 0.00 0.00 0.00 0.00 0.00 0.00 [2,] 0 0.15 0.00 0.00 0.00 0.00 0.00 [3,] 0 0.00 0.15 0.00 0.00 0.00 0.00 [4,] 0 0.00 0.00 0.15 0.00 0.00 0.00 [5,] 0 0.00 0.00 0.00 0.15 0.00 0.00 [6,] 0 0.00 0.00 0.00 0.00 0.15 0.00 [7,] 0 0.00 0.00 0.00 0.00 0.00 0.15 ``` ```r round(unname(cov2cor(vcov(MOP6))), 2) #corr matrix ``` ``` [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1 0 0 0 0 0 0 [2,] 0 1 0 0 0 0 0 [3,] 0 0 1 0 0 0 0 [4,] 0 0 0 1 0 0 0 [5,] 0 0 0 0 1 0 0 [6,] 0 0 0 0 0 1 0 [7,] 0 0 0 0 0 0 1 ``` ]] ]] --- layout: false class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ ## Lung data: compare polynomial regression ]] .row[.content[ ```r print(summary(MP4),4) ``` ``` Call: lm(formula = FEV ~ Age + Age2 + Age3 + Age4, data = dat) Residuals: Min 1Q Median 3Q Max -1.20203 -0.27466 0.02218 0.24525 0.93709 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.5685906 1.0434415 3.420 0.000709 *** Age -1.3941151 0.4529027 -3.078 0.002267 ** Age2 0.2712825 0.0692350 3.918 0.000110 *** Age3 -0.0181507 0.0044318 -4.096 5.37e-05 *** Age4 0.0004055 0.0001007 4.029 7.05e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3897 on 313 degrees of freedom Multiple R-squared: 0.6404, Adjusted R-squared: 0.6358 F-statistic: 139.3 on 4 and 313 DF, p-value: < 2.2e-16 ``` ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ ## Lung data: compare orthogonal polynomial regression ]] .row[.split-70[ .column[.content[ ```r MOP4 <- lm(FEV ~ poly(Age, degree = 4), data = dat) print(summary(MOP4),4) ``` ``` Call: lm(formula = FEV ~ poly(Age, degree = 4), data = dat) Residuals: Min 1Q Median 3Q Max -1.20203 -0.27466 0.02218 0.24525 0.93709 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.45117 0.02185 112.159 < 2e-16 *** poly(Age, degree = 4)1 8.49647 0.38972 21.802 < 2e-16 *** poly(Age, degree = 4)2 -3.14027 0.38972 -8.058 1.65e-14 *** poly(Age, degree = 4)3 -0.35524 0.38972 -0.912 0.363 poly(Age, degree = 4)4 1.56999 0.38972 4.029 7.05e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.3897 on 313 degrees of freedom Multiple R-squared: 0.6404, Adjusted R-squared: 0.6358 F-statistic: 139.3 on 4 and 313 DF, p-value: < 2.2e-16 ``` ]] .column.bg-brand-gray[.content[ * Same model fit from `\(R^2\)` and F-test statistic between MP4 and MOP4. * Note that the p-values does not change much between MOP6 and MOP4 as each order is independent now. * We do not give higher order `\(P_k(x)\)` (at least 2) to avoid complication. * Hence the interpretation of orthogonal polynomial is not straightforward after transformation. ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Response surface model ]] .row[.split-two[ .column[.content[ * Polynomial regression can be extended to models with .brand-blue[two or more covariates]. * The resulting models are often referred to as .brand-blue[response surface models]. * For example, suppose that we have explanatory variables `\(x_{\cdot 1}\)` and `\(x_{\cdot 2}\)` (not data point). * Then a .brand-blue[quadratic response surface model] is `$$E[Y_i] = \beta_0 + \beta_1 x_{i1} + \beta_2 x_ {i2} + \beta_{3} x_{i1} x_{i2} + \beta_4 x_{i1}^2 + \beta_5 x_{i2}^2.$$` ]] .column.bg-brand-gray[.content[ ### Nitrogen oxide emissions for an engine * Ethanol fuel was burned in a single-cylinder engine. * The emission of nitrogen oxide (response variable, `NOx`) was recorded for various settings of the explanatory variables engine compression, `C`, and equivalence ratio, `E`. * Aim in this example is to develop a quadratic response model for `NOx` as a function of `C` and `E`. .bottom_abs.width100[ Reference: Brinkman, N.D. (1981) Ethanol Fuel: A Single-Cylinder Engine Study of Efficiency and Exhaust Emissions. *SAE transactions*, **90**, 1410-1424. ] ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Nitrogen oxide emissions for an engine ]] .row[.split-60[ .column[.content[ ```r dat <- read.table(file = "data/engine.txt", header=TRUE) head(dat,4) ``` ``` NOx C E 1 3.741 12 0.907 2 2.295 12 0.761 3 1.498 12 1.108 4 2.881 12 1.016 ``` ```r summary(dat) ``` ``` NOx C E Min. :0.370 Min. : 7.500 Min. :0.5350 1st Qu.:0.953 1st Qu.: 8.625 1st Qu.:0.7618 Median :1.754 Median :12.000 Median :0.9320 Mean :1.957 Mean :12.034 Mean :0.9265 3rd Qu.:3.003 3rd Qu.:15.000 3rd Qu.:1.1098 Max. :4.028 Max. :18.000 Max. :1.2320 ``` ]] .column.bg-brand-gray[.content[ ![](lecture12_2020JC_files/figure-html/unnamed-chunk-26-1.png)<!-- --> ]] ]] --- class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Engine: Quadratic surface model ]] .row[.split-60[ .column[.content[ ```r M1=lm(NOx ~ C + E + I(C^2) + I(E^2) + I(C * E),data=dat) summary(M1) ``` ``` Call: lm(formula = NOx ~ C + E + I(C^2) + I(E^2) + I(C * E), data = dat) Residuals: Min 1Q Median 3Q Max -0.83859 -0.39735 -0.02915 0.40881 0.76624 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -24.262183 1.545126 -15.702 < 2e-16 *** C 0.224139 0.119449 1.876 0.064152 . E 56.757050 2.756243 20.592 < 2e-16 *** I(C^2) 0.002500 0.004426 0.565 0.573701 I(E^2) -29.785011 1.387610 -21.465 < 2e-16 *** I(C * E) -0.238886 0.061757 -3.868 0.000219 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.4505 on 82 degrees of freedom Multiple R-squared: 0.8509, Adjusted R-squared: 0.8418 F-statistic: 93.61 on 5 and 82 DF, p-value: < 2.2e-16 ``` ]] .column.bg-brand-gray[.content[ ![](lecture12_2020JC_files/figure-html/unnamed-chunk-28-1.png)<!-- --> * One may drop the higher order term .brand-blue[C^2] but not .brand-blue[C]. ]] ]] --- layout: false class: split-10 .row.bg-brand-blue.white[.content.vmiddle[ # Comments on response surface model ]] .row[.content[ * This model has `\(p=2\)` `\(x\)` variables at `\(k=2\)` order. * For `\(p=3\)` and `\(k=2\)`, the full model contains terms `\(x_1\)`, `\(x_2\)`, `\(x_3\)`, `\(x_1^2\)`, `\(x_2^2\)`, `\(x_2^2\)`, `\(x_1x_2\)`, `\(x_1 x_3\)`, `\(x_2 x_3\)`. * The number of terms increases rapidly with increasing `\(p\)` and `\(k\)` and they are highly correlated and becomes a machine learning problem. ]]