# Linear Regression Workflow

This example shows how to fit a linear regression model. A typical workflow involves the following: import data, fit a regression, test its quality, modify it to improve the quality, and share it.

### Step 1. Import the data into a table.

`hospital.xls` is an Excel® spreadsheet containing patient names, sex, age, weight, blood pressure, and dates of treatment in an experimental protocol. First read the data into a table.

`patients = readtable('hospital.xls','ReadRowNames',true);`

Examine the five rows of data.

`patients(1:5,:)`
```ans=5×11 table name sex age wgt smoke sys dia trial1 trial2 trial3 trial4 ____________ _____ ___ ___ _____ ___ ___ ______ ______ ______ ______ YPL-320 {'SMITH' } {'m'} 38 176 1 124 93 18 -99 -99 -99 GLI-532 {'JOHNSON' } {'m'} 43 163 0 109 77 11 13 22 -99 PNI-258 {'WILLIAMS'} {'f'} 38 131 0 125 83 -99 -99 -99 -99 MIJ-579 {'JONES' } {'f'} 40 133 0 117 75 6 12 -99 -99 XLK-030 {'BROWN' } {'f'} 49 119 0 122 80 14 23 -99 -99 ```

The `sex` and `smoke` fields seem to have two choices each. So change these fields to categorical.

```patients.smoke = categorical(patients.smoke,0:1,{'No','Yes'}); patients.sex = categorical(patients.sex);```

### Step 2. Create a fitted model.

Your goal is to model the systolic pressure as a function of a patient's age, weight, sex, and smoking status. Create a linear formula for `'sys'` as a function of `'age'`, `'wgt'`, `'sex'`, and `'smoke'` .

```modelspec = 'sys ~ age + wgt + sex + smoke'; mdl = fitlm(patients,modelspec)```
```mdl = Linear regression model: sys ~ 1 + sex + age + wgt + smoke Estimated Coefficients: Estimate SE tStat pValue _________ ________ ________ __________ (Intercept) 118.28 7.6291 15.504 9.1557e-28 sex_m 0.88162 2.9473 0.29913 0.76549 age 0.08602 0.06731 1.278 0.20438 wgt -0.016685 0.055714 -0.29947 0.76524 smoke_Yes 9.884 1.0406 9.498 1.9546e-15 Number of observations: 100, Error degrees of freedom: 95 Root Mean Squared Error: 4.81 R-squared: 0.508, Adjusted R-Squared: 0.487 F-statistic vs. constant model: 24.5, p-value = 5.99e-14 ```

The sex, age, and weight predictors have rather high $p$-values, indicating that some of these predictors might be unnecessary.

### Step 3. Locate and remove outliers.

See if there are outliers in the data that should be excluded from the fit. Plot the residuals.

`plotResiduals(mdl)` There is one possible outlier, with a value greater than 12. This is probably not truly an outlier. For demonstration, here is how to find and remove it.

Find the outlier.

```outlier = mdl.Residuals.Raw > 12; find(outlier)```
```ans = 84 ```

Remove the outlier.

```mdl = fitlm(patients,modelspec,... 'Exclude',84); mdl.ObservationInfo(84,:)```
```ans=1×4 table Weights Excluded Missing Subset _______ ________ _______ ______ WXM-486 1 true false false ```

Observation 84 is no longer in the model.

### Step 4. Simplify the model.

Try to obtain a simpler model, one with fewer predictors but the same predictive accuracy. `step` looks for a better model by adding or removing one term at a time. Allow `step` take up to 10 steps.

`mdl1 = step(mdl,'NSteps',10)`
```1. Removing wgt, FStat = 4.6001e-05, pValue = 0.9946 2. Removing sex, FStat = 0.063241, pValue = 0.80199 ```
```mdl1 = Linear regression model: sys ~ 1 + age + smoke Estimated Coefficients: Estimate SE tStat pValue ________ ________ ______ __________ (Intercept) 115.11 2.5364 45.383 1.1407e-66 age 0.10782 0.064844 1.6628 0.09962 smoke_Yes 10.054 0.97696 10.291 3.5276e-17 Number of observations: 99, Error degrees of freedom: 96 Root Mean Squared Error: 4.61 R-squared: 0.536, Adjusted R-Squared: 0.526 F-statistic vs. constant model: 55.4, p-value = 1.02e-16 ```

`step` took two steps. This means it could not improve the model further by adding or subtracting a single term.

Plot the effectiveness of the simpler model on the training data.

`plotResiduals(mdl1)` The residuals look about as small as those of the original model.

### Step 5. Predict responses to new data.

Suppose you have four new people, aged 25, 30, 40, and 65, and the first and third smoke. Predict their systolic pressure using `mdl1`.

```ages = [25;30;40;65]; smoker = {'Yes';'No';'Yes';'No'}; systolicnew = feval(mdl1,ages,smoker)```
```systolicnew = 4×1 127.8561 118.3412 129.4734 122.1149 ```

To make predictions, you need only the variables that `mdl1` uses.

### Step 6. Share the model.

You might want others to be able to use your model for prediction. Access the terms in the linear model.

`coefnames = mdl1.CoefficientNames`
```coefnames = 1x3 cell {'(Intercept)'} {'age'} {'smoke_Yes'} ```

View the model formula.

`mdl1.Formula`
```ans = sys ~ 1 + age + smoke ```

Access the coefficients of the terms.

`coefvals = mdl1.Coefficients(:,1).Estimate`
```coefvals = 3×1 115.1066 0.1078 10.0540 ```

The model is `sys = 115.1066 + 0.1078*age + 10.0540*smoke`, where `smoke` is `1` for a smoker, and `0` otherwise.