• R project
• RStudio

• Base R
• Data import
• Data transformation with dplyr
• Data tidying
• Data visualization with ggplot2
• RMarkdown
• Reticulate
• All Cheatsheets

url="https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv"
wine = read.csv2(url, dec = ".")
wine %>% summary

 fixed.acidity   volatile.acidity  citric.acid    residual.sugar  
 Min.   : 4.60   Min.   :0.1200   Min.   :0.000   Min.   : 0.900  
 1st Qu.: 7.10   1st Qu.:0.3900   1st Qu.:0.090   1st Qu.: 1.900  
 Median : 7.90   Median :0.5200   Median :0.260   Median : 2.200  
 Mean   : 8.32   Mean   :0.5278   Mean   :0.271   Mean   : 2.539  
 3rd Qu.: 9.20   3rd Qu.:0.6400   3rd Qu.:0.420   3rd Qu.: 2.600  
 Max.   :15.90   Max.   :1.5800   Max.   :1.000   Max.   :15.500  
   chlorides       free.sulfur.dioxide total.sulfur.dioxide    density      
 Min.   :0.01200   Min.   : 1.00       Min.   :  6.00       Min.   :0.9901  
 1st Qu.:0.07000   1st Qu.: 7.00       1st Qu.: 22.00       1st Qu.:0.9956  
 Median :0.07900   Median :14.00       Median : 38.00       Median :0.9968  
 Mean   :0.08747   Mean   :15.87       Mean   : 46.47       Mean   :0.9967  
 3rd Qu.:0.09000   3rd Qu.:21.00       3rd Qu.: 62.00       3rd Qu.:0.9978  
 Max.   :0.61100   Max.   :72.00       Max.   :289.00       Max.   :1.0037  
       pH          sulphates         alcohol         quality     
 Min.   :2.740   Min.   :0.3300   Min.   : 8.40   Min.   :3.000  
 1st Qu.:3.210   1st Qu.:0.5500   1st Qu.: 9.50   1st Qu.:5.000  
 Median :3.310   Median :0.6200   Median :10.20   Median :6.000  
 Mean   :3.311   Mean   :0.6581   Mean   :10.42   Mean   :5.636  
 3rd Qu.:3.400   3rd Qu.:0.7300   3rd Qu.:11.10   3rd Qu.:6.000  
 Max.   :4.010   Max.   :2.0000   Max.   :14.90   Max.   :8.000  

R has 4 functions for each probability distribution:

• Functions starting with d: probability density
• Functions starting with p: cumulative probability
• Functions starting with r: probability generation
• Functions starting with q: probability quantiles

• Two parameters:
  • Minimum value
  • Maximum value

u = runif(1000)
par(mfrow=c(1,2))
hist(u, prob=T)
curve(dunif(x,0,1),add=T, col="red")
curve(punif(x), col="blue")

summary(u)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0001104 0.2454357 0.4973468 0.4989982 0.7537833 0.9992668 

• Continuous variable
• Most frequently used distribution in statistics and data analysis
• Parameters:
  • \(\mu\) – Mean
  • \(\sigma\) – Standard deviation

x = rnorm(100)
mean(x)

[1] -0.0671774

sd(x)

[1] 1.046495

y = rnorm(100, mean=10, sd = 3)
mean(y)

[1] 9.845252

sd(y)

[1] 3.071808

pnorm(1)

[1] 0.8413447

qnorm(0.975)

[1] 1.959964

hist(x,probability=T,col=gray(.9), ylim=c(0,0.45))
curve(dnorm(x),-5,5, add=T, col="red")

x.all = seq(-3, 3, by = 0.01); y.all = dnorm(x.all)
plot(x.all, y.all, lwd = 2, type = "l", col="red")
x.p = seq(-3, qnorm(.05), length = 100); y.p = dnorm(x.p)
polygon(c(-3, x.p, qnorm(.05)), c(0, y.p, 0), col = "gray")
x.p.up = seq(qnorm(.95), 3, length = 100); y.p.up = dnorm(x.p.up)
polygon(c(qnorm(.95), x.p.up, 3), c(0, y.p.up, 0), col = "gray")

• Continuous variable
• Similar to normal distribution but with longer tails
• Parameter: number of degrees of freedom v
• Approximates normal distribution for large values of v

curve(dnorm, -5, 5, col="blue")
curve(dt(x, df=30), -5, 5, col="purple", add=T)
curve(dt(x, df=10), -5, 5, col="red", add=T)
curve(dt(x, df=1), -5, 5, col="orange", add=T)

• When a sample size, with known mean and variance, grows
• The sample distribution approximates a normal distribution with:
  • Mean equal to the mean of the original distribution
  • Standard deviation equal to the standard deviation of the original distribution divided by the square root of the sample size

\(\mu = \mu_{sample}\)

\(\sigma = \frac{\sigma_{sample}}{\sqrt{n}}\)

Show that the distribution of the mean of samples of size n from a normal distribution with \(\mu = 0\) and \(\sigma = 1\) has standard deviation

\(\sigma = \frac{1}{\sqrt{n}}\)

numsims = 10000; tam = 10; values = rnorm(numsims * tam)
m = matrix(values, numsims, tam)
cat(mean(m), sd(m))

-0.002291752 1.000896

medias = apply(m, 1, mean)
cat(mean(medias), sd(medias))

-0.002291752 0.3169819

1 / sqrt(tam)

[1] 0.3162278

Show that the distribution of the mean of samples of size n from a uniform distribution \(U[0, 1]\) has standard deviation

\(\sigma = \frac{1}{\sqrt{12 \times n}}\)

values.u = runif(numsims * tam)
m.u = matrix(values.u, numsims, tam)
cat(mean(m.u), sd(m.u))

0.4987695 0.2885937

1/sqrt(12)

[1] 0.2886751

medias.u = apply(m.u, 1, mean)
cat(mean(medias.u), sd(medias.u))

0.4987695 0.09188324

1/sqrt(12*tam)

[1] 0.09128709

• In statistical analysis, it is often important to know whether a dataset follows the normal distribution
• An alternative to histograms is to use quantile-to-quantile plots (QQ plots)
• For comparing two distributions, R has the qqplot function
• In order to represent the quantiles of a dataset vs the theoretical quantiles of the normal distribution, we may use the functions qqnorm and qqline

set.seed(123); par(mfrow=c(1,2))
x = rnorm(100, 10, 5); y = runif(100)
qqnorm(x,main='N(10,5)', pch=19, col="blue")
qqline(x, col="red", lwd=3)
qqnorm(y,main='U(0,1)', pch=19, col="blue")
qqline(y, col="red", lwd=3)

• One of the main tasks in data analysis is to estimate the parameters of a distribution from data
• E.g., estimate the mean of a distribution, assuming it is normal
• We cannot find the parameter value but we may estimate it, and quantify the uncertainty
• Confidence intervals provide an estimate, given a confidence level \(\alpha\)

\(P(\hat{\theta_I} < \theta < \hat{\theta_S}) = 1 - \alpha\)

• Estimate the mean of a set of values, assuming they follow the normal distribution
• Usually, we don’t know the standard deviation, thus it will be estimated from the sample
• We use the T statistic that:
  • For small sample size n, follows the t-student distribution with \(n-1\) degrees of freedom
  • For large sample sizes, it follows the normal distribution

\(t = \frac{\bar{X} - \mu}{\frac{s}{\sqrt{n}}}\)

where:

• \(\bar{X}\) is the sample mean
• \(s\) is the sample standard deviation

With known sd

\(\bar{x} - z_{1-\alpha/2} \frac{\sigma}{\sqrt{n}} < \mu < \bar{x} + z_{1-\alpha/2} \frac{\sigma}{\sqrt{n}}\)

With unknown sd and small \(n\)

\(P(\bar{x} - t_{\alpha/2,n-1} \frac{s}{\sqrt{n}} < \mu < \bar{x} + t_{\alpha/2,n-1} \frac{s}{\sqrt{n}}) = 1 - \alpha\)

With unknown sd and large \(n\)

\(\bar{x} - z_{1-\alpha/2} \frac{s}{\sqrt{n}} < \mu < \bar{x} + z_{1-\alpha/2} \frac{s}{\sqrt{n}}\)

x = c(175,176,173,175,174,173,173,176,173,179)
mean(x)

[1] 174.7

se = sd(x) / sqrt(10)
se

[1] 0.6155395

qt(0.975,9)

[1] 2.262157

mean(x) + c(-se*qt(0.975,9), se*qt(0.975,9))

[1] 173.3076 176.0924

t.test(x)

    One Sample t-test

data:  x
t = 283.82, df = 9, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 173.3076 176.0924
sample estimates:
mean of x 
    174.7 

• We formulate a null hypothesis about the parameter values
• \(p\) value is the probability of obtaining the estimate if we assume the null hypothesis
• Sufficiently small values of \(p\) allow us to reject the null hypothesis
• If we don’t reject the null hypothesis, the test is inconclusive

Parametric tests

• They assume the normal distribution
• They focus on specified parameters

Non-parametric tests

• They have no assumption about the data distribution
• There are no parameters

• Null hypothesis: mean is zero (by default)
• Function: t.test
• alternative: two.sided, greater, less

t.test(x, mu=175)

    One Sample t-test

data:  x
t = -0.48738, df = 9, p-value = 0.6376
alternative hypothesis: true mean is not equal to 175
95 percent confidence interval:
 173.3076 176.0924
sample estimates:
mean of x 
    174.7 

wilcox.test(x, conf.int=T)

Warning in wilcox.test.default(x, conf.int = T): cannot compute exact p-value
with ties

Warning in wilcox.test.default(x, conf.int = T): cannot compute exact confidence
interval with ties

    Wilcoxon signed rank test with continuity correction

data:  x
V = 55, p-value = 0.005541
alternative hypothesis: true location is not equal to 0
95 percent confidence interval:
 173 176
sample estimates:
(pseudo)median 
         174.5 

successes = rbinom(1, size = 100, prob = .4)
prop.test(successes, 100, p = 0.6)

    1-sample proportions test with continuity correction

data:  successes out of 100, null probability 0.6
X-squared = 7.5938, df = 1, p-value = 0.005857
alternative hypothesis: true p is not equal to 0.6
95 percent confidence interval:
 0.3608720 0.5622361
sample estimates:
   p 
0.46 

• Null hypothesis is usually the same estimate in both populations
• We use the same functions

• H0: mean(drug) – mean(placebo) = 0
• H1: mean(drug) – mean(placebo) < 0
• Assume equal variances

drug = c(15, 10, 13, 7, 9, 8, 21, 9, 14, 8)
placebo = c(15, 14, 12, 8, 14, 7, 16, 10, 15, 12)
t.test(drug, placebo, alt="less", var.equal=TRUE)

    Two Sample t-test

data:  drug and placebo
t = -0.53311, df = 18, p-value = 0.3002
alternative hypothesis: true difference in means is less than 0
95 percent confidence interval:
     -Inf 2.027436
sample estimates:
mean of x mean of y 
     11.4      12.3 

• In this case, we don’t reject H0

• H0: mean(diag_med1) – mean(diag_med2) = 0
• H1: mean(diag_med1) – mean(diag_med2) != 0
• Assume paired samples (each value of one sample corresponds to the corresponding element in the other sample)

diag_med1 = c(3, 0, 5, 2, 5, 5, 5, 4, 4, 5)
diag_med2 = c(2, 1, 4, 1, 4, 3, 3, 2, 3, 5)
t.test(diag_med1, diag_med2, paired=T)

    Paired t-test

data:  diag_med1 and diag_med2
t = 3.3541, df = 9, p-value = 0.008468
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.325555 1.674445
sample estimates:
mean of the differences 
                      1 

• Reject H0

drug = c(15, 10, 13, 7, 9, 8, 21, 9, 14, 8)
placebo = c(15, 14, 12, 8, 14, 7, 16, 10, 15, 12)
wilcox.test(drug, placebo)

Warning in wilcox.test.default(drug, placebo): cannot compute exact p-value with
ties

    Wilcoxon rank sum test with continuity correction

data:  drug and placebo
W = 39, p-value = 0.4246
alternative hypothesis: true location shift is not equal to 0

• H0: The null hypothesis is that the four populations from which the patients were drawn have the same true proportion of smokers.
• H1: The alternative is that this proportion is different in at least one of the populations.

smokers  <- c( 83, 90, 129, 70 )
patients <- c( 86, 93, 136, 82 )
prop.test(smokers, patients)

    4-sample test for equality of proportions without continuity
    correction

data:  smokers out of patients
X-squared = 12.6, df = 3, p-value = 0.005585
alternative hypothesis: two.sided
sample estimates:
   prop 1    prop 2    prop 3    prop 4 
0.9651163 0.9677419 0.9485294 0.8536585 

• Check whether data comes from a given population
• Applied to categorical variables
• Null hypothesis: each category has probability \(p_i\)

die_freq = c(22,21,22,27,22,36)
probs = rep(1/6,6)
chisq.test(die_freq, p=probs)

    Chi-squared test for given probabilities

data:  die_freq
X-squared = 6.72, df = 5, p-value = 0.2423

letter_freqs = c(100,110,80,55,14)
probs = c(29, 21, 17, 17, 16)/100
chisq.test(letter_freqs, p=probs)

    Chi-squared test for given probabilities

data:  letter_freqs
X-squared = 55.395, df = 4, p-value = 2.685e-11

set.seed(123); sz = 50
A = sample.int(5, sz, replace = TRUE)
B = ifelse(runif(sz) < 0.4, A, sample.int(5, sz, replace = TRUE))
cont.tab = table(A, B)
cont.tab

   B
A   1 2 3 4 5
  1 6 2 2 2 0
  2 1 5 0 0 3
  3 0 2 8 2 0
  4 0 0 1 4 1
  5 0 2 2 2 5

summary(cont.tab)

Number of cases in table: 50 
Number of factors: 2 
Test for independence of all factors:
    Chisq = 50.96, df = 16, p-value = 1.611e-05
    Chi-squared approximation may be incorrect

• Non-parametric test
• Check whether a sample follows a given theoretical distribution
• Check whether two samples come from the same distribution

x=rnorm(1000)
ks.test(x, "pnorm")

    One-sample Kolmogorov-Smirnov test

data:  x
D = 0.023979, p-value = 0.6133
alternative hypothesis: two-sided

ks.test(as.vector(scale(iris$Petal.Length)), "pnorm")

Warning in ks.test(as.vector(scale(iris$Petal.Length)), "pnorm"): ties should
not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  as.vector(scale(iris$Petal.Length))
D = 0.19815, p-value = 1.532e-05
alternative hypothesis: two-sided

ks.test(as.vector(scale(iris$Sepal.Length)), "pnorm")

Warning in ks.test(as.vector(scale(iris$Sepal.Length)), "pnorm"): ties should
not be present for the Kolmogorov-Smirnov test

    One-sample Kolmogorov-Smirnov test

data:  as.vector(scale(iris$Sepal.Length))
D = 0.088654, p-value = 0.1891
alternative hypothesis: two-sided

shapiro.test(rnorm(100))

    Shapiro-Wilk normality test

data:  rnorm(100)
W = 0.99156, p-value = 0.7884

shapiro.test(runif(100))

    Shapiro-Wilk normality test

data:  runif(100)
W = 0.95099, p-value = 0.0009647

par(mfrow=c(1,2))
hist(iris$Sepal.Width, col = "cyan", breaks = 10, prob = T)
lines(density(iris$Sepal.Width), col = "red", lwd = 2)
qqnorm(iris$Sepal.Width)
qqline(iris$Sepal.Width, col = "red", lwd = 2)

shapiro.test(iris$Sepal.Width)

    Shapiro-Wilk normality test

data:  iris$Sepal.Width
W = 0.98492, p-value = 0.1012

• If a test has \(\alpha = 1\%\), the expected number of type I errors is 1 in 100
• By doing \(N\) tests for a given \(\alpha\), we will have of average \((1-\alpha) \times N\) false positives!
• We need to perform \(p\) value corrections!

• Control FWER (family wise error rate)
  • Define a new \(\alpha\) by dividing the original value by the number of tests (Bonferroni correction)
  • Very conservative method
• Control FDR (false discovery rate)
  • corrections BH (Benjamini / Hochberg) e BY
  • Less conservative
• Function p.adjust makes the p-value correction according to one of these methods

set.seed(12345)
m = matrix(rnorm(100*60),100,60)
tt = function(x) t.test(x[1:30],x[31:60])$p.value
pvalues = apply(m, 1, tt)
sum(pvalues < 0.05)

[1] 5

adj.bonf = p.adjust(pvalues, "bonferroni")
sum(adj.bonf < 0.05)

[1] 0

adj.bh = p.adjust(pvalues, "BH")
sum(adj.bh < 0.05)

[1] 0

set.seed(12345)
# m1: all come from population with the same mean (50 instances)
m1 = matrix(rnorm(50*60),50,60)
# m2: half with mean 0 and half with mean 1
m2 = cbind(matrix(rnorm(50*30),50,30), matrix(rnorm(50*30, mean = 1), 50, 30))
mt = rbind(m1, m2)
pvalues.2 = apply(mt, 1, tt)
res.verd = c(rep("don't rej",50), rep("rej",50))
table(pvalues.2 < 0.05, res.verd)

       res.verd
        don't rej rej
  FALSE        47   3
  TRUE          3  47

adj.bonf.2 = p.adjust(pvalues.2, "bonferroni")
table(adj.bonf.2 < 0.05, res.verd)

       res.verd
        don't rej rej
  FALSE        50  20
  TRUE          0  30

adj.bh.2 = p.adjust(pvalues.2, "BH")
table(adj.bh.2 < 0.05, res.verd)

       res.verd
        don't rej rej
  FALSE        49   4
  TRUE          1  46

• Analysis of variance is a set of statistical models and tests that aims to partition the observational variance in several components
• Simpler version aims to understand whether the mean of several groups is the same
• Generalizes tests to the means for \(N > 2\) samples
• Statistical test that assumes that data follows the normal distribution with null hypothesis that the means are the same in all groups

prof1 = c(4,3,4,5,2,3,4,5); prof2 = c(4,4,5,5,4,5,4,4); prof3 = c(3,4,2,4,5,5,4,4)
notas = data.frame(prof1, prof2, prof3)
notas_st = stack(notas)
notas_st %>% group_by(ind) %>% summarize(
  Count = length(values),
  Mean = mean(values),
  SD = sd(values))

# A tibble: 3 x 4
  ind   Count  Mean    SD
  <fct> <int> <dbl> <dbl>
1 prof1     8  3.75 1.04 
2 prof2     8  4.38 0.518
3 prof3     8  3.88 0.991

1. Perform a test to check if variances are different
2. Perform a test to check for differences in groups

fligner.test(values ~ ind, data= notas_st)

    Fligner-Killeen test of homogeneity of variances

data:  values by ind
Fligner-Killeen:med chi-squared = 1.3956, df = 2, p-value = 0.4977

Because fligner.test didn’t indicate that variances were different we indicate that var.equal is true

oneway.test(values ~ ind, data=notas_st, var.equal = TRUE)

    One-way analysis of means

data:  values and ind
F = 1.1308, num df = 2, denom df = 21, p-value = 0.3417

• Kruskal-Wallis Rank Sum Test
• Non-parametric test
• Does not depend on the normality of data
• Assumes same shaped distributions
• Null hypothesis is that all groups have the same median
• The alternative is that at least one group is different

set.seed(42)
df = stack(data.frame(p1 = rnorm(20, mean = 20, sd = 4),
         p2 = rnorm(20, mean = 25, sd = 3), p3 = rnorm(20, mean = 27, sd = 3)))

kruskal.test(values ~ ind, data=df)

    Kruskal-Wallis rank sum test

data:  values by ind
Kruskal-Wallis chi-squared = 17.965, df = 2, p-value = 0.0001256

oneway.test(values ~ ind, data=df)

    One-way analysis of means (not assuming equal variances)

data:  values and ind
F = 11.113, num df = 2.00, denom df = 36.64, p-value = 0.0001689

TukeyHSD(aov(values ~ ind, data=df))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = values ~ ind, data = df)

$ind
          diff        lwr      upr     p adj
p2-p1 3.419345  0.3729782 6.465711 0.0242937
p3-p1 6.229588  3.1832214 9.275954 0.0000227
p3-p2 2.810243 -0.2361231 5.856610 0.0763708

plot(TukeyHSD(aov(values ~ ind, data=df)))

pairwise.t.test(df$values, df$ind, p.adjust.method = "bonferroni")

    Pairwise comparisons using t tests with pooled SD 

data:  df$values and df$ind 

   p1      p2   
p2 0.027   -    
p3 2.3e-05 0.091

P value adjustment method: bonferroni 

pairwise.wilcox.test(df$values, df$ind, p.adjust.method = "bonferroni")

    Pairwise comparisons using Wilcoxon rank sum exact test 

data:  df$values and df$ind 

   p1      p2     
p2 0.05861 -      
p3 0.00027 0.01166

P value adjustment method: bonferroni 

head(mtcars, 2)

              mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4      21   6  160 110  3.9 2.620 16.46  0  1    4    4
Mazda RX4 Wag  21   6  160 110  3.9 2.875 17.02  0  1    4    4

round(cor(mtcars), 2)

       mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
mpg   1.00 -0.85 -0.85 -0.78  0.68 -0.87  0.42  0.66  0.60  0.48 -0.55
cyl  -0.85  1.00  0.90  0.83 -0.70  0.78 -0.59 -0.81 -0.52 -0.49  0.53
disp -0.85  0.90  1.00  0.79 -0.71  0.89 -0.43 -0.71 -0.59 -0.56  0.39
hp   -0.78  0.83  0.79  1.00 -0.45  0.66 -0.71 -0.72 -0.24 -0.13  0.75
drat  0.68 -0.70 -0.71 -0.45  1.00 -0.71  0.09  0.44  0.71  0.70 -0.09
wt   -0.87  0.78  0.89  0.66 -0.71  1.00 -0.17 -0.55 -0.69 -0.58  0.43
qsec  0.42 -0.59 -0.43 -0.71  0.09 -0.17  1.00  0.74 -0.23 -0.21 -0.66
vs    0.66 -0.81 -0.71 -0.72  0.44 -0.55  0.74  1.00  0.17  0.21 -0.57
am    0.60 -0.52 -0.59 -0.24  0.71 -0.69 -0.23  0.17  1.00  0.79  0.06
gear  0.48 -0.49 -0.56 -0.13  0.70 -0.58 -0.21  0.21  0.79  1.00  0.27
carb -0.55  0.53  0.39  0.75 -0.09  0.43 -0.66 -0.57  0.06  0.27  1.00

heatmap(cor(mtcars))

library(corrplot)
corrplot(cor(mtcars),method = "number", type = "upper")

x = c(0.3, 0.8, 0.5, 0.5, 0.2, 0.2)
y = c(0.3, 0.6, 0.5, 0.8, 0.2, 0.1)
cor.test(x, y)

    Pearson's product-moment correlation

data:  x and y
t = 2.4893, df = 4, p-value = 0.06753
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.08711455  0.97455731
sample estimates:
      cor 
0.7795607 

#library(Hmisc)
res = rcorr(as.matrix(mtcars)) # rcorr() accepts matrices only
round(res$P, 3)

       mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
mpg     NA 0.000 0.000 0.000 0.000 0.000 0.017 0.000 0.000 0.005 0.001
cyl  0.000    NA 0.000 0.000 0.000 0.000 0.000 0.000 0.002 0.004 0.002
disp 0.000 0.000    NA 0.000 0.000 0.000 0.013 0.000 0.000 0.001 0.025
hp   0.000 0.000 0.000    NA 0.010 0.000 0.000 0.000 0.180 0.493 0.000
drat 0.000 0.000 0.000 0.010    NA 0.000 0.620 0.012 0.000 0.000 0.621
wt   0.000 0.000 0.000 0.000 0.000    NA 0.339 0.001 0.000 0.000 0.015
qsec 0.017 0.000 0.013 0.000 0.620 0.339    NA 0.000 0.206 0.243 0.000
vs   0.000 0.000 0.000 0.000 0.012 0.001 0.000    NA 0.357 0.258 0.001
am   0.000 0.002 0.000 0.180 0.000 0.000 0.206 0.357    NA 0.000 0.754
gear 0.005 0.004 0.001 0.493 0.000 0.000 0.243 0.258 0.000    NA 0.129
carb 0.001 0.002 0.025 0.000 0.621 0.015 0.000 0.001 0.754 0.129    NA

#library(GGally)
ggpairs(mtcars[, c("mpg", "hp", "wt", "cyl")])

set.seed(1234)
x = rnorm(100)
e = rnorm(100, mean=0, sd=5)
y = 2 - 7 * x + e
plot(x, y, pch = 16, cex = 1.3, col = "blue")
abline(lm(y ~ x))

lm(y ~ x) %>% summary

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.4313  -3.0700   0.0118   2.9322  14.9387 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.1858     0.5249   4.164 6.73e-05 ***
x            -7.1304     0.5189 -13.742  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.186 on 98 degrees of freedom
Multiple R-squared:  0.6583,    Adjusted R-squared:  0.6549 
F-statistic: 188.8 on 1 and 98 DF,  p-value: < 2.2e-16

lm(y ~ x) %>% confint

                2.5 %    97.5 %
(Intercept)  1.144136  3.227405
x           -8.160124 -6.100725

• Residuals should follow a normal distribution

par(mfrow=c(2,2))
lm(y ~ x) %>% plot

• Points 78, 91 and 92 are probably outliers

lm(y ~ x) %>% outlierTest

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
   rstudent unadjusted p-value Bonferroni p
78 3.017626          0.0032541      0.32541

lm(y ~ x) %>% influence.measures

Influence measures of
     lm(formula = y ~ x) :

       dfb.1_     dfb.x     dffit cov.r   cook.d    hat inf
1    2.80e-02 -3.56e-02  0.049198 1.040 1.22e-03 0.0210    
2   -5.18e-02 -2.13e-02 -0.053481 1.028 1.44e-03 0.0119    
3    6.64e-03  6.98e-03  0.008966 1.047 4.06e-05 0.0254    
4   -3.97e-02  1.34e-01 -0.147585 1.076 1.10e-02 0.0580   *
5   -8.97e-02 -4.87e-02 -0.096364 1.020 4.66e-03 0.0134    
6    1.52e-02  9.23e-03  0.016707 1.035 1.41e-04 0.0144    
7   -8.54e-02  3.87e-02 -0.100223 1.015 5.03e-03 0.0117    
8    1.05e-02 -4.42e-03  0.012165 1.032 7.47e-05 0.0115    
9    2.72e-02 -1.20e-02  0.031786 1.031 5.10e-04 0.0117    
10  -9.58e-03  8.04e-03 -0.013592 1.036 9.33e-05 0.0154    
11  -2.24e-02  7.64e-03 -0.025018 1.031 3.16e-04 0.0110    
12  -5.97e-02  5.87e-02 -0.091106 1.028 4.17e-03 0.0171    
13  -1.02e-01  7.08e-02 -0.134465 1.008 9.01e-03 0.0138    
14   8.09e-02  1.75e-02  0.081062 1.018 3.30e-03 0.0105    
15   1.17e-03  1.12e-03  0.001506 1.044 1.15e-06 0.0225    
16   7.65e-02  3.58e-03  0.076973 1.019 2.97e-03 0.0100    
17  -1.18e-01  4.49e-02 -0.134364 0.999 8.97e-03 0.0113    
18   9.18e-03 -7.95e-03  0.013201 1.037 8.80e-05 0.0157    
19   5.29e-02 -4.08e-02  0.072466 1.028 2.64e-03 0.0146    
20  -5.96e-05 -1.11e-04 -0.000119 1.105 7.10e-09 0.0763   *
21  -2.26e-02 -6.35e-03 -0.022736 1.031 2.61e-04 0.0108    
22  -7.58e-02  2.71e-02 -0.085366 1.019 3.66e-03 0.0111    
23   1.88e-01 -5.64e-02  0.206622 0.953 2.07e-02 0.0108    
24   7.66e-02  4.36e-02  0.083104 1.025 3.47e-03 0.0138    
25   1.58e-01 -9.38e-02  0.198155 0.973 1.92e-02 0.0129    
26   3.98e-04 -6.53e-04  0.000826 1.049 3.45e-07 0.0267    
27  -7.03e-02 -4.67e-02 -0.079073 1.028 3.15e-03 0.0154    
28  -1.33e-01  1.35e-01 -0.206349 0.990 2.10e-02 0.0175    
29  -6.61e-02 -9.27e-03 -0.066085 1.022 2.20e-03 0.0102    
30   1.39e-02 -1.25e-02  0.020368 1.037 2.10e-04 0.0161    
31   1.18e-01  1.25e-01  0.160026 1.027 1.28e-02 0.0259    
32   1.85e-02 -6.29e-03  0.020692 1.031 2.16e-04 0.0110    
33  -1.08e-01  6.65e-02 -0.137351 1.004 9.39e-03 0.0131    
34   5.62e-02 -2.07e-02  0.063584 1.025 2.03e-03 0.0112    
35  -1.32e-01  2.56e-01 -0.309706 0.993 4.70e-02 0.0317    
36  -3.53e-02  4.30e-02 -0.060421 1.038 1.84e-03 0.0202    
37  -2.80e-02  8.40e-02 -0.093693 1.072 4.43e-03 0.0510   *
38  -1.63e-01  2.40e-01 -0.314214 0.965 4.79e-02 0.0240    
39   8.23e-02 -1.17e-02  0.085939 1.016 3.70e-03 0.0102    
40  -6.14e-02  2.02e-02 -0.068392 1.023 2.35e-03 0.0110    
41  -4.11e-02 -5.33e-02 -0.062823 1.056 1.99e-03 0.0358    
42   1.11e-01 -1.20e-01  0.177480 1.005 1.56e-02 0.0183    
43   4.95e-02 -3.94e-02  0.068690 1.030 2.38e-03 0.0149    
44  -1.44e-02  1.84e-03 -0.014934 1.031 1.13e-04 0.0102    
45   3.78e-02 -3.69e-02  0.057442 1.034 1.66e-03 0.0170    
46   2.83e-02 -2.67e-02  0.042357 1.036 9.05e-04 0.0166    
47   1.33e-01 -1.50e-01  0.218219 0.990 2.35e-02 0.0190    
48   1.65e-02 -2.22e-02  0.029996 1.043 4.54e-04 0.0220    
49   4.12e-02 -1.62e-02  0.047126 1.028 1.12e-03 0.0113    
50   2.70e-02 -9.82e-03  0.030493 1.030 4.69e-04 0.0112    
51  -3.37e-02  7.60e-02 -0.088825 1.056 3.98e-03 0.0372    
52   4.05e-03 -1.87e-03  0.004777 1.033 1.15e-05 0.0118    
53   1.31e-01 -1.48e-01  0.215012 0.991 2.28e-02 0.0191    
54  -7.87e-02  7.91e-02 -0.121378 1.021 7.38e-03 0.0174    
55   7.69e-03 -4.32e-05  0.007788 1.031 3.06e-05 0.0100    
56   1.45e-01  9.48e-02  0.162229 1.001 1.31e-02 0.0152    
57  -2.91e-02 -4.14e-02 -0.047331 1.065 1.13e-03 0.0426   *
58  -9.70e-02  6.71e-02 -0.127710 1.011 8.14e-03 0.0138    
59  -1.10e-01 -1.53e-01 -0.176118 1.049 1.56e-02 0.0411    
60  -3.73e-02  4.49e-02 -0.063450 1.037 2.03e-03 0.0200    
61  -9.46e-02 -6.91e-02 -0.109492 1.023 6.01e-03 0.0166    
62  -3.41e-02 -6.56e-02 -0.069881 1.112 2.47e-03 0.0833   *
63  -4.42e-02 -5.36e-03 -0.044220 1.027 9.86e-04 0.0101    
64  -2.10e-02  1.19e-02 -0.025962 1.033 3.40e-04 0.0126    
65   3.63e-02  5.35e-03  0.036252 1.029 6.63e-04 0.0102    
66   8.24e-02  1.24e-01  0.139316 1.063 9.77e-03 0.0474   *
67   1.33e-01 -1.57e-01  0.223868 0.990 2.47e-02 0.0197    
68  -8.52e-03 -1.06e-02 -0.012701 1.056 8.15e-05 0.0333    
69  -3.91e-02 -4.77e-02 -0.057459 1.052 1.67e-03 0.0321    
70   1.51e-01  6.99e-02  0.157856 0.993 1.23e-02 0.0124    
71   1.66e-01  2.68e-02  0.165986 0.977 1.35e-02 0.0103    
72  -5.29e-04  1.68e-04 -0.000586 1.032 1.74e-07 0.0109    
73  -3.52e-02  7.72e-03 -0.037604 1.029 7.13e-04 0.0104    
74  -2.03e-01 -1.47e-01 -0.234640 0.971 2.69e-02 0.0165    
75   1.96e-01  3.28e-01  0.359567 1.041 6.40e-02 0.0597    
76  -8.45e-02 -2.88e-04 -0.085465 1.016 3.66e-03 0.0100    
77  -9.45e-02  1.46e-01 -0.188521 1.018 1.77e-02 0.0252    
78   2.73e-01 -1.72e-01  0.349236 0.864 5.63e-02 0.0132   *
79   2.09e-02  8.25e-03  0.021509 1.032 2.34e-04 0.0117    
80  -7.34e-03  1.22e-03 -0.007724 1.031 3.01e-05 0.0103    
81  -2.75e-01  5.88e-03 -0.279354 0.885 3.65e-02 0.0100   *
82  -1.35e-02  1.81e-04 -0.013702 1.031 9.48e-05 0.0100    
83   7.13e-02 -1.08e-01  0.140437 1.030 9.88e-03 0.0248    
84   3.56e-02 -6.15e-04  0.036093 1.028 6.57e-04 0.0100    
85   1.01e-01  8.89e-02  0.125205 1.026 7.86e-03 0.0202    
86   2.19e-01  1.67e-01  0.257173 0.963 3.22e-02 0.0173    
87   1.23e-01  7.95e-02  0.137608 1.010 9.44e-03 0.0150    
88  -5.13e-02  1.33e-02 -0.055639 1.026 1.56e-03 0.0106    
89   6.32e-02 -2.24e-03  0.064324 1.022 2.08e-03 0.0100    
90  -2.16e-02  2.71e-02 -0.037683 1.041 7.17e-04 0.0208    
91  -5.62e-02 -5.80e-03 -0.056268 1.025 1.59e-03 0.0101    
92  -3.07e-01 -1.20e-01 -0.315916 0.874 4.64e-02 0.0117   *
93  -1.01e-01 -1.47e-01 -0.166679 1.056 1.39e-02 0.0447    
94   5.47e-02  5.43e-02  0.071773 1.041 2.60e-03 0.0234    
95   1.96e-01 -7.12e-02  0.221776 0.945 2.38e-02 0.0111    
96   4.91e-02  2.36e-02  0.051677 1.029 1.35e-03 0.0126    
97   4.54e-02 -5.31e-02  0.075876 1.035 2.90e-03 0.0196    
98  -1.11e-01 -9.99e-02 -0.138872 1.023 9.65e-03 0.0207    
99   1.72e-02  1.67e-02  0.022356 1.044 2.52e-04 0.0228    
100 -2.89e-01 -4.91e-01 -0.535996 0.997 1.39e-01 0.0620   *

set.seed(42)
u = rnorm(100)
v = rnorm(100, mean = 3, sd = 2)
w = rnorm(100, mean = -3, sd = 1)
e = rnorm(100, mean = 0, sd = 3)
y = 5 + 4 * u + 3 * v + 2 * w + e
df = data.frame(u, v, w, y)
lm(y ~ u + v + w, data = df)

Call:
lm(formula = y ~ u + v + w, data = df)

Coefficients:
(Intercept)            u            v            w  
      4.770        4.173        3.013        1.905  

lm(y ~ u + v + w, data = df) %>% summary

Call:
lm(formula = y ~ u + v + w, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3832 -1.7601 -0.3115  1.8565  6.9840 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.7703     0.9691   4.922 3.55e-06 ***
u             4.1733     0.2597  16.069  < 2e-16 ***
v             3.0132     0.1484  20.311  < 2e-16 ***
w             1.9052     0.2665   7.150 1.71e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.66 on 96 degrees of freedom
Multiple R-squared:  0.8851,    Adjusted R-squared:  0.8815 
F-statistic: 246.6 on 3 and 96 DF,  p-value: < 2.2e-16

• The residuals should have a perfect normal distribution with mean of zero
• The sign of the median indicates the skew’s direction
• The magnitude of the median indicates the extent
• There is a skew to the left
• The 1st and 3rd quantiles should have about the same magnitude
• In this case, there is a slight skew to the right
• The min and max residuals allow us to detect extreme outliers

lm(y ~ u + v + w, data = df) %>% residuals %>% summary

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-5.3832 -1.7601 -0.3115  0.0000  1.8565  6.9840 

lm(y ~ u + v + w, data = df) %>% anova

Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value    Pr(>F)    
u          1 1776.24 1776.24 251.005 < 2.2e-16 ***
v          1 3097.08 3097.08 437.655 < 2.2e-16 ***
w          1  361.74  361.74  51.118 1.707e-10 ***
Residuals 96  679.35    7.08                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Add interactions between factors
• Use multiplication instead of additions in the formula
• \(u * v\) is the same as \(u + v + u:v\)
• \(u * v * w\) is the same as \(u + v + u:v + u:w + v:w + u:v:w\)
• (u + v + w)^2 is the same as \(u + v + u:v + u:w + v:w\)

set.seed(42)
u = rnorm(100); v = rnorm(100, mean = 3, sd = 2); w = rnorm(100, mean = -2, sd = 1)
e = rnorm(100, mean = 0, sd = 5)
y = 5 + 4 * u - 3 * v * w + e

lm(y ~ u + v + w) %>% summary

Call:
lm(formula = y ~ u + v + w)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.7515  -3.4697  -0.3039   4.4165  14.0817 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -11.8346     1.8321  -6.460 4.32e-09 ***
u             5.3375     0.6323   8.441 3.28e-13 ***
v             5.1933     0.3612  14.378  < 2e-16 ***
w            -9.4179     0.6488 -14.516  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.477 on 96 degrees of freedom
Multiple R-squared:  0.842, Adjusted R-squared:  0.837 
F-statistic: 170.5 on 3 and 96 DF,  p-value: < 2.2e-16

lm(y ~ u * v * w) %>% summary

Call:
lm(formula = y ~ u * v * w)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.1595 -3.0684 -0.3264  3.0443 11.8704 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.0941     2.1213   1.930   0.0567 .  
u             3.6410     1.8798   1.937   0.0558 .  
v             0.2456     0.5606   0.438   0.6623    
w            -0.5976     1.0169  -0.588   0.5582    
u:v           0.1302     0.3908   0.333   0.7398    
u:w          -0.5478     0.8700  -0.630   0.5305    
v:w          -2.8597     0.2832 -10.099   <2e-16 ***
u:v:w         0.1540     0.2188   0.704   0.4834    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.496 on 92 degrees of freedom
Multiple R-squared:  0.927, Adjusted R-squared:  0.9215 
F-statistic: 166.9 on 7 and 92 DF,  p-value: < 2.2e-16

step(lm(y ~ 1), direction="forward", scope = ( ~ u * v * w), trace = 0) %>% summary

Call:
lm(formula = y ~ w + v + u + w:v)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0948 -2.9481 -0.3681  3.1480 11.5484 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.5983     1.9396   1.855   0.0667 .  
w            -0.8142     0.9362  -0.870   0.3867    
v             0.3885     0.5225   0.744   0.4590    
u             4.3631     0.4436   9.836  3.7e-16 ***
w:v          -2.7874     0.2669 -10.445  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.442 on 95 degrees of freedom
Multiple R-squared:  0.9264,    Adjusted R-squared:  0.9233 
F-statistic: 299.1 on 4 and 95 DF,  p-value: < 2.2e-16

step(lm(y ~ (u + v + w)^2), direction="backward", trace = 0) %>% summary

Call:
lm(formula = y ~ u + v + w + v:w)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0948 -2.9481 -0.3681  3.1480 11.5484 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.5983     1.9396   1.855   0.0667 .  
u             4.3631     0.4436   9.836  3.7e-16 ***
v             0.3885     0.5225   0.744   0.4590    
w            -0.8142     0.9362  -0.870   0.3867    
v:w          -2.7874     0.2669 -10.445  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.442 on 95 degrees of freedom
Multiple R-squared:  0.9264,    Adjusted R-squared:  0.9233 
F-statistic: 299.1 on 4 and 95 DF,  p-value: < 2.2e-16

step(lm(y ~ u * v * w), direction="backward", trace = 0) %>% summary

Call:
lm(formula = y ~ u + v + w + v:w)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0948 -2.9481 -0.3681  3.1480 11.5484 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.5983     1.9396   1.855   0.0667 .  
u             4.3631     0.4436   9.836  3.7e-16 ***
v             0.3885     0.5225   0.744   0.4590    
w            -0.8142     0.9362  -0.870   0.3867    
v:w          -2.7874     0.2669 -10.445  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.442 on 95 degrees of freedom
Multiple R-squared:  0.9264,    Adjusted R-squared:  0.9233 
F-statistic: 299.1 on 4 and 95 DF,  p-value: < 2.2e-16

set.seed(4)
n <- 150
x1 <- rnorm(n)
x2 <- rnorm(n, 1, 2)
x3 <- rnorm(n, 3, 1)
x4 <- rnorm(n,-2, 2)
e <- rnorm(n, 0, 3)
y <- 4 + x1 + 5 * x3 + e

step(lm(y ~ x1*x2*x3*x4), direction="backward", trace = 0) %>% summary

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x1:x2 + x1:x3 + x2:x3 + 
    x1:x4 + x2:x4 + x1:x2:x3 + x1:x2:x4)

Residuals:
   Min     1Q Median     3Q    Max 
-7.933 -1.770  0.157  2.220  6.823 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.27344    0.88242   3.710   0.0003 ***
x1           0.13632    0.94325   0.145   0.8853    
x2          -0.04601    0.43371  -0.106   0.9157    
x3           5.18404    0.27093  19.134   <2e-16 ***
x4          -0.03134    0.15519  -0.202   0.8402    
x1:x2       -0.22227    0.45474  -0.489   0.6258    
x1:x3        0.11347    0.28704   0.395   0.6932    
x2:x3       -0.02440    0.13040  -0.187   0.8519    
x1:x4       -0.01450    0.14950  -0.097   0.9229    
x2:x4       -0.06272    0.07980  -0.786   0.4332    
x1:x2:x3     0.27541    0.14192   1.941   0.0543 .  
x1:x2:x4     0.13491    0.08075   1.671   0.0971 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.893 on 138 degrees of freedom
Multiple R-squared:  0.7854,    Adjusted R-squared:  0.7682 
F-statistic:  45.9 on 11 and 138 DF,  p-value: < 2.2e-16

step(lm(y ~ 1), direction="forward", scope = (~ x1*x2*x3*x4), trace = 0) %>% summary

Call:
lm(formula = y ~ x3 + x1)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.1478 -1.8497 -0.0554  2.0257  6.5499 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.6481     0.7510   4.857 3.01e-06 ***
x3            5.1469     0.2386  21.568  < 2e-16 ***
x1            0.5824     0.2552   2.282   0.0239 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.923 on 147 degrees of freedom
Multiple R-squared:  0.7666,    Adjusted R-squared:  0.7635 
F-statistic: 241.5 on 2 and 147 DF,  p-value: < 2.2e-16

library(purrr)
df = data.frame(x1,x2,x3,x4,y)
df %>% select(-y) %>% map_dbl(cor, y = df$y) %>% sort(decreasing = TRUE)

          x3           x1           x2           x4 
 0.870840575  0.167670468  0.001738779 -0.021349925 

df %>% select(-y) %>% map_dbl(cor, y = df$y) %>% sort(decreasing = TRUE) %>% .[1:2]

       x3        x1 
0.8708406 0.1676705 

df %>% select(-y) %>% map_dbl(cor, y = df$y) %>% sort(decreasing = TRUE) %>% .[1:2] %>% names 

[1] "x3" "x1"

best2 = df %>% select(-y) %>% map_dbl(cor, y = df$y) %>% sort(dec = TRUE) %>% .[1:2] %>% names %>% df[.]
lm(df$y ~ as.matrix(best2)) %>% summary

Call:
lm(formula = df$y ~ as.matrix(best2))

Residuals:
    Min      1Q  Median      3Q     Max 
-8.1478 -1.8497 -0.0554  2.0257  6.5499 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)          3.6481     0.7510   4.857 3.01e-06 ***
as.matrix(best2)x3   5.1469     0.2386  21.568  < 2e-16 ***
as.matrix(best2)x1   0.5824     0.2552   2.282   0.0239 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.923 on 147 degrees of freedom
Multiple R-squared:  0.7666,    Adjusted R-squared:  0.7635 
F-statistic: 241.5 on 2 and 147 DF,  p-value: < 2.2e-16

model = lm(y ~ x1 + x3, data = df)
vals = data.frame(x1 = c(1, 0.8), x2 = c(1, 1.2), x3 = c(4, 3.7), x4 = c(-1.8,-2))
predict(model, newdata = vals)

       1        2 
24.81809 23.15755 

predict(model, newdata = vals, interval="prediction")

       fit      lwr      upr
1 24.81809 18.98696 30.64923
2 23.15755 17.34215 28.97295