This code examines the association of five genes with age (KRT14, CCL21, CD1a, CD3E, CXCL12) adjusted by KRT8. We used univariate analysis of variance (ANOVA) to evaluate the association.
The associated file “GeneExpressionFileCreation.rmd” replaces expression values that are below the limit of detection (LOD) for each gene with 1/2 the LOD. It also cleans up names and id’s so that it can be merged with the subject characteristics file (e.g., age and percent thymic epithelium levels)
require(dplyr)
require(plyr)
require(rms)
require(readr)
raw.data <- read_csv("Quantgene2.csv")
## data formats
final.data <- as.data.frame(raw.data) %>%
mutate(age.ch=as.factor(ifelse(`Age of Thymus Tissue Donor (in Years)`<=17, "Young","Adult"))) %>% # categorizing age to two groups in age.ch variable
mutate(age.ch3=as.factor(ifelse(`Age of Thymus Tissue Donor (in Years)`<=1, "Postnatal",
ifelse(`Age of Thymus Tissue Donor (in Years)`>1 & `Age of Thymus Tissue Donor (in Years)` <=17,"Child-Adolescent","Adult")))) %>%
# categorizing age to two or three groups in age.ch and age.ch3 variable
rename(replace = c("Total area of relevant tissue on slide:" = "relevant.tissue")) %>%
rename(replace = c("Total area of lymphoid tissue on slide:" = "lymphoid.tissue")) %>%
rename(replace = c("Total area of thymic epithelium on slide:" = "thymic.epithelium")) %>%
rename(replace = c("Total active cortical area (thymopoiesis):" = "cortical")) %>%
rename(replace=c("DHVI Thymus Tissue ID Number (T#)"="id")) %>%
rename(replace =c("Age of Thymus Tissue Donor (in Years)"="Age"))
names(final.data) <- gsub(" \\(human\\)","",names(final.data))
names(final.data) <- gsub(" ","",names(final.data))
#create \% thymic epithelial and \% cortical
area.to.total <- final.data %>%
mutate(lymphoid.tissue.p=100*lymphoid.tissue/relevant.tissue,
thymic.epithelium.p=100*thymic.epithelium/relevant.tissue,
cortical.p=100*cortical/relevant.tissue) %>%
select(id, Age,age.ch,age.ch3,Gender,
lymphoid.tissue.p, thymic.epithelium.p, cortical.p,
everything(),
-relevant.tissue, -lymphoid.tissue,-thymic.epithelium, -cortical)
area.to.total$age.ch3 <- factor(area.to.total$age.ch3, levels=c("Postnatal", "Child-Adolescent", "Adult"))
new.expression <-read.csv("new.expression.csv", stringsAsFactors=F)
### P10 = P1, 9086= 908 missing 244 -- can I assume 244 was a typo in the original clinical #####
sort(new.expression$id)
## [1] "1013" "1080" "158" "2000" "2001" "2002" "2003" "222" "240" "263"
## [11] "264" "265" "278" "309" "321" "349" "359" "545" "562" "577"
## [21] "587" "588" "591" "598" "628" "630" "639" "680" "697" "707"
## [31] "717" "908" "921" "924" "926" "960" "997" "P1" "P11" "P12"
## [41] "P13" "P14" "P15" "P16" "P17" "P18" "P19" "P2" "P20" "P21"
## [51] "P22" "P23" "P24" "P25" "P4" "P5" "P6" "P7" "P8" "P9"
sort(area.to.total$id)
## [1] "1013" "1080" "158" "2000" "2001" "2002" "2003" "222" "240" "244"
## [11] "265" "278" "309" "321" "349" "359" "545" "562" "577" "587"
## [21] "588" "591" "598" "628" "630" "639" "680" "697" "707" "717"
## [31] "908" "921" "924" "926" "960" "997" "P1" "P11" "P13" "P17"
## [41] "P2" "P20" "P21" "P22" "P24" "P4" "P8"
names(area.to.total)
## [1] "id" "Age" "age.ch"
## [4] "age.ch3" "Gender" "lymphoid.tissue.p"
## [7] "thymic.epithelium.p" "cortical.p" "VCAM1"
## [10] "IL12a" "PPARG" "KRT14"
## [13] "LIF" "AIRE" "NCAM1"
## [16] "CCL25" "PSMB11" "FLT3LG"
## [19] "FGF7" "BMP4" "IL6"
## [22] "ADIPOQ" "KRT8" "IL22"
## [25] "CD1a" "GAPDH" "IL23a"
## [28] "KIT" "CCL19" "ICAM1"
## [31] "DLL4" "IL7" "FOXP3"
## [34] "ATF3" "S1PR1" "ESR1"
## [37] "CCL21" "AR" "IN3C1"
## [40] "GUSB" "HPRT1" "RORC"
## [43] "TNFSF11" "CD3E" "WNTF4"
## [46] "ITGAX" "KITLG" "NOTCH1"
## [49] "IL2" "FOXN1" "CCND1"
## [52] "CXCL12" "GREM1" "CD80"
## [55] "JAG1" "TNFRSF11b" "LTB"
## [58] "ITH3" "CYP21a2" "SELP"
dim(area.to.total)
## [1] 47 60
area.to.total2 <-subset(area.to.total, select = 1:8)
updated.data <-merge(area.to.total2, new.expression, by="id")
updated.data$agenew <-ifelse(updated.data$Age >= 60,"60 +", updated.data$age.ch)
knitr::kable(table(updated.data$agenew))
Var1 | Freq |
---|---|
1 | 16 |
2 | 25 |
60 + | 5 |
knitr::kable(table(updated.data$Gender))
Var1 | Freq |
---|---|
Female | 27 |
Male | 18 |
Unknown | 1 |
### Show that GAPDH expression changes with age
GAPDH <- lm(log(GAPDH) ~ age.ch, data=updated.data)
summary(GAPDH)
##
## Call:
## lm(formula = log(GAPDH) ~ age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8009 -0.7738 -0.4675 0.4242 3.9292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.3965 0.3172 17.011 <2e-16 ***
## age.chYoung 0.9160 0.4303 2.129 0.0389 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.454 on 44 degrees of freedom
## Multiple R-squared: 0.09337, Adjusted R-squared: 0.07276
## F-statistic: 4.531 on 1 and 44 DF, p-value: 0.03892
GAPDH.mean <- lm(log(GAPDH) ~ -1 + age.ch, data=updated.data)
summary(GAPDH.mean)
##
## Call:
## lm(formula = log(GAPDH) ~ -1 + age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8009 -0.7738 -0.4675 0.4242 3.9292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age.chAdult 5.3965 0.3172 17.01 <2e-16 ***
## age.chYoung 6.3124 0.2907 21.71 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.454 on 44 degrees of freedom
## Multiple R-squared: 0.9453, Adjusted R-squared: 0.9428
## F-statistic: 380.4 on 2 and 44 DF, p-value: < 2.2e-16
### Does KRT14 show similar difference in age?
KRT14 <- lm(log(KRT14) ~ age.ch, data=updated.data)
summary(KRT14)
##
## Call:
## lm(formula = log(KRT14) ~ age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9581 -0.8388 0.0583 0.9165 3.4747
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.2162 0.2927 21.236 <2e-16 ***
## age.chYoung 0.4154 0.3971 1.046 0.301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.341 on 44 degrees of freedom
## Multiple R-squared: 0.02427, Adjusted R-squared: 0.002094
## F-statistic: 1.094 on 1 and 44 DF, p-value: 0.3012
KRT14means <-lm(log(KRT14) ~ -1 + age.ch, data=updated.data)
summary(KRT14means)
##
## Call:
## lm(formula = log(KRT14) ~ -1 + age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9581 -0.8388 0.0583 0.9165 3.4747
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age.chAdult 6.2162 0.2927 21.24 <2e-16 ***
## age.chYoung 6.6316 0.2683 24.72 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.341 on 44 degrees of freedom
## Multiple R-squared: 0.9602, Adjusted R-squared: 0.9584
## F-statistic: 531 on 2 and 44 DF, p-value: < 2.2e-16
##### Does KRT8 show similar difference in age?
KRT8<- lm(log(KRT8) ~ age.ch, data=updated.data)
summary(KRT8)
##
## Call:
## lm(formula = log(KRT8) ~ age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8360 -0.5644 -0.2303 0.2739 3.0758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0665 0.1879 26.97 <2e-16 ***
## age.chYoung 0.2879 0.2548 1.13 0.265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8608 on 44 degrees of freedom
## Multiple R-squared: 0.0282, Adjusted R-squared: 0.006114
## F-statistic: 1.277 on 1 and 44 DF, p-value: 0.2646
KRT8means <-lm(log(KRT8) ~ -1 + age.ch, data=updated.data)
summary(KRT8means)
##
## Call:
## lm(formula = log(KRT8) ~ -1 + age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8360 -0.5644 -0.2303 0.2739 3.0758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age.chAdult 5.0665 0.1879 26.97 <2e-16 ***
## age.chYoung 5.3545 0.1722 31.10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8608 on 44 degrees of freedom
## Multiple R-squared: 0.9747, Adjusted R-squared: 0.9735
## F-statistic: 847.3 on 2 and 44 DF, p-value: < 2.2e-16
########################
KRT8 <- lm(log(KRT8/KRT14) ~ age.ch, data=updated.data)
anova(KRT8)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 0.1854245 | 0.1854245 | 0.1528252 | 0.6977365 |
Residuals | 44 | 53.3856894 | 1.2133111 | NA | NA |
####################
CD1a <- lm(log(CD1a/KRT14) ~ age.ch, data=updated.data)
anova(CD1a)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 2.798422 | 2.798422 | 2.19558 | 0.1455339 |
Residuals | 44 | 56.081109 | 1.274571 | NA | NA |
####################
CCL21 <- lm(log(CCL.21/ KRT14) ~ age.ch, data=updated.data)
anova(CCL21)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 4.702175 | 4.702175 | 2.481426 | 0.1223615 |
Residuals | 44 | 83.377747 | 1.894949 | NA | NA |
####################
CD3E <- lm(log(CD3.E/KRT14) ~ age.ch, data=updated.data)
anova(CD3E)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 18.53158 | 18.5315797 | 29.25367 | 2.474859e-06 |
Residuals | 44 | 27.87307 | 0.6334788 | NA | NA |
####################
CXCL12 <- lm(log(CXCL.12/KRT14) ~ age.ch, data=updated.data)
anova(CXCL12)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 0.929426 | 0.9294260 | 0.998508 | 0.3231357 |
Residuals | 44 | 40.955853 | 0.9308148 | NA | NA |
########################
KRT14 <- lm(log(KRT14/KRT8) ~ age.ch, data=updated.data)
anova(KRT14)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 0.1854245 | 0.1854245 | 0.1528252 | 0.6977365 |
Residuals | 44 | 53.3856894 | 1.2133111 | NA | NA |
KRT14.means <- lm(KRT14/KRT8 ~ -1 + age.ch, data=updated.data)
summary(KRT14.means)
##
## Call:
## lm(formula = KRT14/KRT8 ~ -1 + age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.643 -2.928 -1.246 1.486 13.529
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age.chAdult 5.8017 0.9686 5.990 3.49e-07 ***
## age.chYoung 4.7378 0.8878 5.337 3.15e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.439 on 44 degrees of freedom
## Multiple R-squared: 0.5939, Adjusted R-squared: 0.5755
## F-statistic: 32.18 on 2 and 44 DF, p-value: 2.449e-09
####################
CD1a.krt8 <- lm(log(CD1a/KRT8) ~ age.ch, data=updated.data)
anova(CD1a.krt8)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 4.424535 | 4.4245349 | 5.241784 | 0.02690493 |
Residuals | 44 | 37.139942 | 0.8440896 | NA | NA |
CD1a.krt8.means <- lm(CD1a/KRT8 ~ -1 + age.ch, data=updated.data)
summary(CD1a.krt8.means)
##
## Call:
## lm(formula = CD1a/KRT8 ~ -1 + age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9927 -1.2898 -0.8416 0.8922 7.9366
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age.chAdult 1.7287 0.5108 3.384 0.00151 **
## age.chYoung 3.1636 0.4682 6.758 2.59e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.341 on 44 degrees of freedom
## Multiple R-squared: 0.5649, Adjusted R-squared: 0.5451
## F-statistic: 28.56 on 2 and 44 DF, p-value: 1.121e-08
####################
CCL21.krt8 <- lm(log(CCL.21/ KRT8) ~ age.ch, data=updated.data)
anova(CCL21.krt8)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 3.020089 | 3.020089 | 3.03736 | 0.08835384 |
Residuals | 44 | 43.749815 | 0.994314 | NA | NA |
CCL21.krt8.means <- lm(CCL.21/ KRT8 ~ -1 + age.ch, data=updated.data)
summary(CCL21.krt8.means)
##
## Call:
## lm(formula = CCL.21/KRT8 ~ -1 + age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.984 -2.737 -0.589 0.212 38.194
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age.chAdult 4.264 1.360 3.135 0.00306 **
## age.chYoung 1.128 1.247 0.905 0.37064
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.233 on 44 degrees of freedom
## Multiple R-squared: 0.1948, Adjusted R-squared: 0.1582
## F-statistic: 5.322 on 2 and 44 DF, p-value: 0.008511
####################
CD3E.krt8 <- lm(log(CD3.E/KRT8) ~ age.ch, data=updated.data)
anova(CD3E.krt8)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 22.42441 | 22.424407 | 22.18823 | 2.493195e-05 |
Residuals | 44 | 44.46834 | 1.010644 | NA | NA |
CD3E.krt8.means <- lm(CD3.E/KRT8 ~ -1 + age.ch, data=updated.data)
summary(CD3E.krt8.means)
##
## Call:
## lm(formula = CD3.E/KRT8 ~ -1 + age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2151 -2.2259 -0.9735 1.2725 9.7558
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age.chAdult 2.6502 0.8181 3.240 0.00228 **
## age.chYoung 6.9412 0.7498 9.258 6.82e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.749 on 44 degrees of freedom
## Multiple R-squared: 0.6862, Adjusted R-squared: 0.6719
## F-statistic: 48.1 on 2 and 44 DF, p-value: 8.463e-12
####################
CXCL12.krt8 <- lm(log(CXCL.12/KRT8) ~ age.ch, data=updated.data)
anova(CXCL12.krt8)
Df <int> | Sum Sq <dbl> | Mean Sq <dbl> | F value <dbl> | Pr(>F) <dbl> | |
---|---|---|---|---|---|
age.ch | 1 | 0.2845774 | 0.2845774 | 1.731326 | 0.1950563 |
Residuals | 44 | 7.2322623 | 0.1643696 | NA | NA |
CXCL12.krt8.means <- lm(CXCL.12/KRT8 ~ -1 + age.ch, data=updated.data)
summary(CXCL12.krt8.means)
##
## Call:
## lm(formula = CXCL.12/KRT8 ~ -1 + age.ch, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.986 -2.320 -0.350 2.491 5.035
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## age.chAdult 7.5520 0.5866 12.87 < 2e-16 ***
## age.chYoung 6.2209 0.5376 11.57 6.16e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.688 on 44 degrees of freedom
## Multiple R-squared: 0.872, Adjusted R-squared: 0.8661
## F-statistic: 149.8 on 2 and 44 DF, p-value: < 2.2e-16
library(beanplot)
library(tidyr)
updated.data.nogapdh <-subset(updated.data, select=-c(GAPDH))
updated.data.long <- updated.data.nogapdh %>% gather(Gene, Expression, c(KRT14, CD1a:CXCL.12))
#png("beanplot expression.png")
beanplot(log(Expression/KRT8)~age.ch*Gene, data=updated.data.long, ll = 0.04,
side = "both", xlab="", ylab="Log Ratio to KRT8", cex.axis=0.75,
col = list("grey56", c("grey92", "black")),
axes=T)
#dev.off()
library(dplyr)
library(ggplot2)
library(ggpubr)
theme_set(theme_pubclean())
updated.data.long$Expression.adjusted <-log(updated.data.long$Expression/updated.data.long$KRT8)
#png("boxplot expression.png")
ggplot(updated.data.long, aes(Gene, Expression.adjusted)) +
geom_boxplot(aes(color = age.ch))+ labs(y= "Log Expression Standardized to KRT8") +
scale_color_manual(values = c("black", "grey41")) +
stat_compare_means(aes(group = age.ch), label = "p.signif")
#dev.off()
##
## 0 4.9 7 9 16 17 22 27 28 29 32 34 36 38 39 40 45 49 53 58
## 20 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1
## 59 60 64 67 70 78
## 1 1 1 1 1 1
##
## Call:
## lm(formula = thymic.epithelium.p ~ Gender, data = updated.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.727 -28.767 3.004 30.261 46.019
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.274 6.488 6.670 3.88e-08 ***
## GenderMale 15.897 10.258 1.550 0.129
## GenderUnknown 9.416 34.331 0.274 0.785
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 33.71 on 43 degrees of freedom
## Multiple R-squared: 0.05306, Adjusted R-squared: 0.009021
## F-statistic: 1.205 on 2 and 43 DF, p-value: 0.3097
##
## Call:
## lm(formula = thymic.epithelium.p ~ age.zip, data = area.to.total)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.357 -9.181 3.026 11.506 23.109
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.06707 3.27195 23.25 < 2e-16 ***
## age.zip -1.21148 0.09884 -12.26 6.12e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.49 on 45 degrees of freedom
## Multiple R-squared: 0.7695, Adjusted R-squared: 0.7644
## F-statistic: 150.2 on 1 and 45 DF, p-value: 6.123e-16
##
## Call:
## lm(formula = thymic.epithelium.p ~ Age, data = non.zero)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.896 -9.717 0.108 12.674 34.090
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 61.6117 7.0463 8.744 4.48e-09 ***
## Age -0.9154 0.1613 -5.674 6.59e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.36 on 25 degrees of freedom
## Multiple R-squared: 0.5629, Adjusted R-squared: 0.5454
## F-statistic: 32.19 on 1 and 25 DF, p-value: 6.592e-06
##
## Call:
## lm(formula = thymic.epithelium.p ~ Age, data = area.to.total)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.436 -9.201 2.979 11.495 22.975
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.21779 3.28803 23.18 < 2e-16 ***
## Age -1.21414 0.09932 -12.22 6.72e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.52 on 45 degrees of freedom
## Multiple R-squared: 0.7685, Adjusted R-squared: 0.7634
## F-statistic: 149.4 on 1 and 45 DF, p-value: 6.724e-16
## png
## 2
## png
## 2