4-NP Analysis of %OC
Thurber
November 21, 2018
Import Data
We use the readxl package to import data from sheets in “4NP Standard Curve and Data.xlsx”.
percoc <- read_excel("4NP Standard Curve and Data.xlsx",
sheet = "Data_Complete_%OC",
skip = 2)## New names:
## * `Mass of dust per trial (g)` -> `Mass of dust per trial (g)...5`
## * `4NP in trial (mg)` -> `4NP in trial (mg)...6`
## * `Mass released (mg)` -> `Mass released (mg)...7`
## * `%released` -> `%released...8`
## * `Mass of dust per trial (g)` -> `Mass of dust per trial (g)...9`
## * ... names(percoc)## [1] "Days" "Mass dosed"
## [3] "Mass dust dosed(g)" "Conc on dust (mg/g)"
## [5] "Mass of dust per trial (g)...5" "4NP in trial (mg)...6"
## [7] "Mass released (mg)...7" "%released...8"
## [9] "Mass of dust per trial (g)...9" "4NP in trial (mg)...10"
## [11] "Mass released (mg)...11" "%released...12"
## [13] "Mass of dust per trial (g)...13" "4NP in trial (mg)...14"
## [15] "Mass released (mg)...15" "%released...16"
## [17] "Mass of dust per trial (g)...17" "4NP in trial (mg)...18"
## [19] "Mass released (mg)...19" "%released...20" percoc.4perc <- percoc[, c(1:4,5:8)]
names(percoc.4perc) <- c("Days", "Mass_Dosed", "Dust_Massed_Dosed", "Conc_On_Dust", "Mass_Of_Dust_Per_Trial", "FourNP_In_Trial", "Mass_Released", "Perc_Released")
percoc.4perc$Perc_OC <- 4
percoc.10perc <- percoc[, c(1:4,9:12)]
names(percoc.10perc) <- c("Days", "Mass_Dosed", "Dust_Massed_Dosed", "Conc_On_Dust", "Mass_Of_Dust_Per_Trial", "FourNP_In_Trial", "Mass_Released", "Perc_Released")
percoc.10perc$Perc_OC <- 10
percoc.50perc <- percoc[, c(1:4,13:16)]
names(percoc.50perc) <- c("Days", "Mass_Dosed", "Dust_Massed_Dosed", "Conc_On_Dust", "Mass_Of_Dust_Per_Trial", "FourNP_In_Trial", "Mass_Released", "Perc_Released")
percoc.50perc$Perc_OC <- 50
percoc.70perc <- percoc[, c(1:4,17:20)]
names(percoc.70perc) <- c("Days", "Mass_Dosed", "Dust_Massed_Dosed", "Conc_On_Dust", "Mass_Of_Dust_Per_Trial", "FourNP_In_Trial", "Mass_Released", "Perc_Released")
percoc.70perc$Perc_OC <- 70
percoc <- rbind(percoc.4perc, percoc.10perc, percoc.50perc, percoc.70perc)
percoc$Sample <- factor(round((percoc$Days - floor(percoc$Days))*10))
percoc$Days <- floor(percoc$Days)
percoc$Perc_Released <- as.numeric(percoc$Perc_Released)
percoc$Prop_Released <- percoc$Perc_Released / 100
rm(percoc.4perc, percoc.10perc, percoc.50perc, percoc.70perc)
write.csv(percoc, "Four_NP_PercOC.csv", row.names = FALSE)
head(percoc)## # A tibble: 6 x 11
## Days Mass_Dosed Dust_Massed_Dos~ Conc_On_Dust Mass_Of_Dust_Pe~
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 NA
## 2 5 18.0 5.5 3.27 0.88
## 3 5 15.6 5.6 2.79 1.01
## 4 5 22.3 7 3.19 1.4
## 5 10 18.0 5.5 3.27 1.33
## 6 10 15.6 5.6 2.79 1.21
## # ... with 6 more variables: FourNP_In_Trial <dbl>, Mass_Released <dbl>,
## # Perc_Released <dbl>, Perc_OC <dbl>, Sample <fct>, Prop_Released <dbl> table(percoc$Days, percoc$Perc_OC)##
## 4 10 50 70
## 0 1 1 1 1
## 5 3 3 3 3
## 10 3 3 3 3
## 15 3 3 3 3
## 20 3 3 3 3
## 25 3 3 3 3
## 30 3 3 3 3
## 35 3 3 3 3
## 40 3 3 3 3 ggplot(percoc, aes(x=Days, y=Perc_Released, group=factor(Perc_OC), colour=factor(Perc_OC), shape=factor(Perc_OC), linetype=factor(Perc_OC))) +
geom_point() +
stat_summary(fun.y="mean", geom="line") +
#stat_summary(fun.y="mean", geom="point", shape = "-", size=10) +
theme_bw() +
labs(shape = "OC%", linetype = "OC%", colour = "OC%") +
ylab("Percent 4-NP Released")## Warning: `fun.y` is deprecated. Use `fun` instead.
Fit
Subset the data to get Days = 40.
percoc.40 <- filter(percoc, Days == 40)
write.csv(percoc.40, "Four_NP_PercOC_40.csv", row.names = FALSE)Plot the data to look for trends.
xyplot(Perc_Released ~ Perc_OC, group = Sample, data = percoc.40, ylab = "Percent 4-NP Released", xlab = "Percent OC", cex = 2, pch = 3, col = c("red", "blue", "green"))
Fit a model at Time = 40 using Perc_OC while controlling for Sample.
percoc.40 <- filter(percoc, Days == 40)
percoc.40 %>%
group_by(Perc_OC) %>%
summarize(
avg_perc_oc_released = mean(Perc_Released),
sd_perc_oc_released = sd(Perc_Released)
)## # A tibble: 4 x 3
## Perc_OC avg_perc_oc_released sd_perc_oc_released
## * <dbl> <dbl> <dbl>
## 1 4 13.6 0.231
## 2 10 9.35 0.266
## 3 50 2.1 0.173
## 4 70 0.477 0.0115 ### Unequal sd's indicates possible need for transformation
percoc.40 %>%
group_by(Perc_OC) %>%
summarize(
avg_log_prop_oc_released = mean(log(Prop_Released)),
sd_log_prop_oc_released = sd(log(Prop_Released))
)## # A tibble: 4 x 3
## Perc_OC avg_log_prop_oc_released sd_log_prop_oc_released
## * <dbl> <dbl> <dbl>
## 1 4 -1.99 0.0169
## 2 10 -2.37 0.0287
## 3 50 -3.87 0.0807
## 4 70 -5.35 0.0241 ### Unequal sd's indicates possible need for transformation
percoc.40 %>%
group_by(Perc_OC) %>%
summarize(
avg_sqrt_prop_oc_released = mean(Prop_Released^(1/2)),
sd_sqrt_prop_oc_released = sd(Prop_Released^(1/2))
)## # A tibble: 4 x 3
## Perc_OC avg_sqrt_prop_oc_released sd_sqrt_prop_oc_released
## * <dbl> <dbl> <dbl>
## 1 4 0.369 0.00312
## 2 10 0.306 0.00436
## 3 50 0.145 0.00591
## 4 70 0.0690 0.000833 ### Unequal sd's indicates possible need for transformation
percoc.40$Rank_Prop_Released <- rank(percoc.40$Prop_Released)
percoc.40 %>%
group_by(Perc_OC) %>%
summarize(
avg_log_prop_oc_released = mean(Rank_Prop_Released),
sd_log_prop_oc_released = sd(Rank_Prop_Released)
)## # A tibble: 4 x 3
## Perc_OC avg_log_prop_oc_released sd_log_prop_oc_released
## * <dbl> <dbl> <dbl>
## 1 4 11 0.866
## 2 10 8 0.866
## 3 50 5 0.866
## 4 70 2 0.866 ### Fit the model
percoc.lm <- lm(Perc_Released ~ Perc_OC*Sample, data = percoc.40)
summary(percoc.lm)##
## Call:
## lm(formula = Perc_Released ~ Perc_OC * Sample, data = percoc.40)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.79170 -1.26563 -0.07669 1.13499 1.94509
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.703707 1.518053 8.368 0.000158 ***
## Perc_OC -0.187200 0.035021 -5.345 0.001752 **
## Sample2 -0.013095 2.146851 -0.006 0.995331
## Sample3 -0.013095 2.146851 -0.006 0.995331
## Perc_OC:Sample2 -0.001549 0.049527 -0.031 0.976058
## Perc_OC:Sample3 -0.001549 0.049527 -0.031 0.976058
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.927 on 6 degrees of freedom
## Multiple R-squared: 0.9353, Adjusted R-squared: 0.8813
## F-statistic: 17.33 on 5 and 6 DF, p-value: 0.001653 anova(percoc.lm)## Analysis of Variance Table
##
## Response: Perc_Released
## Df Sum Sq Mean Sq F value Pr(>F)
## Perc_OC 1 321.76 321.76 86.6695 8.696e-05 ***
## Sample 2 0.01 0.01 0.0015 0.9985
## Perc_OC:Sample 2 0.00 0.00 0.0007 0.9993
## Residuals 6 22.27 3.71
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 plot(percoc.lm)
Replot without Sample.
xyplot(Prop_Released ~ Perc_OC, data = percoc.40, ylab = "Percent 4-NP Released")
Since Sample appears to not be a factor, remove it and fit the reduced, complete model.
ggplot(percoc.40, aes(x=Perc_OC, y=Perc_Released)) +
geom_smooth(method = "lm") +
geom_point(shape=1, size=1) +
theme_bw() ## `geom_smooth()` using formula 'y ~ x'
percoc.lm <- lm(Perc_Released ~ Perc_OC, data = percoc.40)
summary(percoc.lm)##
## Call:
## lm(formula = Perc_Released ~ Perc_OC, data = percoc.40)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.77265 -1.29065 -0.01599 1.11800 1.95796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.69498 0.67914 18.69 4.15e-09 ***
## Perc_OC -0.18823 0.01567 -12.01 2.89e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.493 on 10 degrees of freedom
## Multiple R-squared: 0.9352, Adjusted R-squared: 0.9287
## F-statistic: 144.3 on 1 and 10 DF, p-value: 2.889e-07 anova(percoc.lm)## Analysis of Variance Table
##
## Response: Perc_Released
## Df Sum Sq Mean Sq F value Pr(>F)
## Perc_OC 1 321.76 321.76 144.34 2.889e-07 ***
## Residuals 10 22.29 2.23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 plot(percoc.lm)
percoc.lm <- lm(Perc_Released ~ factor(Perc_OC), data = percoc.40)
summary(percoc.lm)##
## Call:
## lm(formula = Perc_Released ~ factor(Perc_OC), data = percoc.40)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.306667 -0.108333 -0.006667 0.153333 0.266667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.6333 0.1133 120.35 2.54e-14 ***
## factor(Perc_OC)10 -4.2867 0.1602 -26.76 4.10e-09 ***
## factor(Perc_OC)50 -11.5333 0.1602 -71.99 1.54e-12 ***
## factor(Perc_OC)70 -13.1567 0.1602 -82.12 5.39e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1962 on 8 degrees of freedom
## Multiple R-squared: 0.9991, Adjusted R-squared: 0.9988
## F-statistic: 2976 on 3 and 8 DF, p-value: 1.58e-12 anova(percoc.lm)## Analysis of Variance Table
##
## Response: Perc_Released
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Perc_OC) 3 343.74 114.579 2976.1 1.58e-12 ***
## Residuals 8 0.31 0.039
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 plot(percoc.lm)
Perc_OC is statistically significant.
The linear model is highly suspect since the scatterplot displays a distinctly nonlinear trend. Interpretation of the ANOVA model is that mean percent released decreases as percent OC increases. However, variability increases as percent released increases (Perc_OC decreases). A transformation may be needed to make variability uniform.
Stupid Jim Assumption: Perc_OC is a continuous covariate and we can get values other than 4, 10, 50, and 70%. So, a transformed model to linearize is preferred over ANOVA. However, with only four observed 0C%s, fitting any function is not going to work very well. It’s probably better to fit the ANOVA and discuss trend over OC%.
We can see the differences in group means by looking at a model on Perc_OC. The Bonferroni, Holm, and Tukey tests and confidence intervals all show the decreasing rate of release as Perc_OC increases.
percoc.aov <- aov(Perc_Released ~ factor(Perc_OC), data = percoc.40)
summary(percoc.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Perc_OC) 3 343.7 114.58 2976 1.58e-12 ***
## Residuals 8 0.3 0.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ### Make room for long y-axis labels
par(mar=c(5,12,2,2)+0.1)
TukeyHSD(percoc.aov)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Perc_Released ~ factor(Perc_OC), data = percoc.40)
##
## $`factor(Perc_OC)`
## diff lwr upr p adj
## 10-4 -4.286667 -4.799710 -3.773624 0.00e+00
## 50-4 -11.533333 -12.046376 -11.020290 0.00e+00
## 70-4 -13.156667 -13.669710 -12.643624 0.00e+00
## 50-10 -7.246667 -7.759710 -6.733624 0.00e+00
## 70-10 -8.870000 -9.383043 -8.356957 0.00e+00
## 70-50 -1.623333 -2.136376 -1.110290 3.59e-05 plot(TukeyHSD(percoc.aov, conf.level=.95), las=1, cex.axis=0.75)
### Reset mar to default
par(mar=c(5,4,4,2)+0.1)
pairwise.t.test(percoc.40$Perc_Released, factor(percoc.40$Perc_OC), p.adj = "bonf")##
## Pairwise comparisons using t tests with pooled SD
##
## data: percoc.40$Perc_Released and factor(percoc.40$Perc_OC)
##
## 4 10 50
## 10 2.5e-08 - -
## 50 9.3e-12 3.8e-10 -
## 70 3.2e-12 7.5e-11 4.6e-05
##
## P value adjustment method: bonferroni pairwise.t.test(percoc.40$Perc_Released, factor(percoc.40$Perc_OC), p.adj = "holm")##
## Pairwise comparisons using t tests with pooled SD
##
## data: percoc.40$Perc_Released and factor(percoc.40$Perc_OC)
##
## 4 10 50
## 10 8.2e-09 - -
## 50 7.7e-12 1.9e-10 -
## 70 3.2e-12 5.0e-11 7.7e-06
##
## P value adjustment method: holmMaybe a logit transformation will stabilize the linear model.
logit <- function (p, percents = range.p[2] > 1, adjust)
{
range.p <- range(p, na.rm = TRUE)
if (percents) {
if (range.p[1] < 0 || range.p[1] > 100)
stop("p must be in the range 0 to 100")
p <- p/100
range.p <- range.p/100
}
else if (range.p[1] < 0 || range.p[1] > 1)
stop("p must be in the range 0 to 1")
a <- if (missing(adjust)) {
if (isTRUE(all.equal(range.p[1], 0)) || isTRUE(all.equal(range.p[2],
1)))
0.025
else 0
}
else adjust
if (missing(adjust) && a != 0)
warning(paste("proportions remapped to (", a, ", ", 1 -
a, ")", sep = ""))
a <- 1 - 2 * a
log((0.5 + a * (p - 0.5))/(1 - (0.5 + a * (p - 0.5))))
}
LogitProp <- logit(percoc.40$Prop_Released)
ggplot(percoc.40, aes(x=Perc_OC, y=LogitProp)) +
geom_smooth(method = "lm") +
geom_point(shape=16, size=1) +
theme_bw() ## `geom_smooth()` using formula 'y ~ x'
percoc.lm <- lm(LogitProp ~ Perc_OC, percoc.40)
summary(percoc.lm)##
## Call:
## lm(formula = LogitProp ~ Perc_OC, data = percoc.40)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21427 -0.15663 -0.04756 0.09023 0.39761
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.660081 0.096744 -17.16 9.54e-09 ***
## Perc_OC -0.049730 0.002232 -22.28 7.44e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2127 on 10 degrees of freedom
## Multiple R-squared: 0.9803, Adjusted R-squared: 0.9783
## F-statistic: 496.5 on 1 and 10 DF, p-value: 7.444e-10 plot(percoc.lm)
While not optimal, we can treat the Perc_OC values as fixed effects and check for homogeneity of variance using Levene’s test.
p_load(car)## Installing package into 'C:/Users/School/OneDrive - University of Redlands/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)## also installing the dependencies 'matrixStats', 'RcppArmadillo', 'numDeriv', 'SparseM', 'MatrixModels', 'conquer', 'minqa', 'nloptr', 'statmod', 'carData', 'abind', 'pbkrtest', 'quantreg', 'maptools', 'rio', 'lme4'## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0:
## cannot open URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0/PACKAGES'##
## There are binary versions available but the source versions are later:
## binary source needs_compilation
## SparseM 1.78 1.81 TRUE
## quantreg 5.83 5.85 TRUE
##
## package 'matrixStats' successfully unpacked and MD5 sums checked
## package 'RcppArmadillo' successfully unpacked and MD5 sums checked
## package 'numDeriv' successfully unpacked and MD5 sums checked
## package 'MatrixModels' successfully unpacked and MD5 sums checked
## package 'conquer' successfully unpacked and MD5 sums checked
## package 'minqa' successfully unpacked and MD5 sums checked
## package 'nloptr' successfully unpacked and MD5 sums checked
## package 'statmod' successfully unpacked and MD5 sums checked
## package 'carData' successfully unpacked and MD5 sums checked
## package 'abind' successfully unpacked and MD5 sums checked
## package 'pbkrtest' successfully unpacked and MD5 sums checked
## package 'maptools' successfully unpacked and MD5 sums checked
## package 'rio' successfully unpacked and MD5 sums checked
## package 'lme4' successfully unpacked and MD5 sums checked
## package 'car' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\School\AppData\Local\Temp\RtmpcP4VSm\downloaded_packages## installing the source packages 'SparseM', 'quantreg'##
## car installed leveneTest(Prop_Released~factor(Perc_OC), data=percoc.40)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.3287 0.805
## 8 leveneTest(logit(Prop_Released)~factor(Perc_OC), data=percoc.40)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.3867 0.7657
## 8According to Levene, whatever lack of homogeneity does exist, it is not of concern.
Box and Cox may provide a nice transformation to normalize the relationship. The EnvStats package is Steve Millard’s and goes with his textbook.
p_load(EnvStats)## Installing package into 'C:/Users/School/OneDrive - University of Redlands/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)## also installing the dependency 'nortest'## Warning: unable to access index for repository http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0:
## cannot open URL 'http://www.stats.ox.ac.uk/pub/RWin/bin/windows/contrib/4.0/PACKAGES'## package 'nortest' successfully unpacked and MD5 sums checked
## package 'EnvStats' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\School\AppData\Local\Temp\RtmpcP4VSm\downloaded_packages##
## EnvStats installed xyplot(Prop_Released ~ Perc_OC, data = percoc.40)
LogitProp <- logit(percoc.40$Prop_Released)
xyplot(LogitProp ~ Perc_OC, data = percoc.40, type="p,r")
summary(percoc.lm <- lm(LogitProp ~ Perc_OC, data = percoc.40))##
## Call:
## lm(formula = LogitProp ~ Perc_OC, data = percoc.40)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21427 -0.15663 -0.04756 0.09023 0.39761
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.660081 0.096744 -17.16 9.54e-09 ***
## Perc_OC -0.049730 0.002232 -22.28 7.44e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2127 on 10 degrees of freedom
## Multiple R-squared: 0.9803, Adjusted R-squared: 0.9783
## F-statistic: 496.5 on 1 and 10 DF, p-value: 7.444e-10 plot(percoc.lm)
### Try Box and Cox
epsilon <- 0.000000001
summary(percoc.lm <- lm((Prop_Released+epsilon) ~ Perc_OC, data = percoc.40))##
## Call:
## lm(formula = (Prop_Released + epsilon) ~ Perc_OC, data = percoc.40)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0177265 -0.0129065 -0.0001599 0.0111800 0.0195796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1269498 0.0067914 18.69 4.15e-09 ***
## Perc_OC -0.0018823 0.0001567 -12.01 2.89e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01493 on 10 degrees of freedom
## Multiple R-squared: 0.9352, Adjusted R-squared: 0.9287
## F-statistic: 144.3 on 1 and 10 DF, p-value: 2.889e-07 boxcox(percoc.lm)## $lambda
## [1] -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0
##
## $objective
## [1] 0.9219219 0.9233033 0.9256661 0.9263268 0.9432943 0.9809693 0.9217730
## [8] 0.9289774 0.9462395
##
## $objective.name
## [1] "PPCC"
##
## $optimize
## [1] FALSE
##
## $optimize.bounds
## lower upper
## NA NA
##
## $eps
## [1] 2.220446e-16
##
## $lm.obj
##
## Call:
## lm(formula = (Prop_Released + epsilon) ~ Perc_OC, data = percoc.40,
## y = TRUE, qr = TRUE)
##
## Coefficients:
## (Intercept) Perc_OC
## 0.126950 -0.001882
##
##
## $sample.size
## [1] 12
##
## $data.name
## [1] "percoc.lm"
##
## attr(,"class")
## [1] "boxcoxLm" box_prop_released <- boxcoxTransform(percoc.40$Prop_Released+epsilon, lambda = 0.5)
summary(percoc.lm <- lm(box_prop_released ~ Perc_OC, data = percoc.40))##
## Call:
## lm(formula = box_prop_released ~ Perc_OC, data = percoc.40)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.047503 -0.021611 0.006348 0.018300 0.044623
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.2641718 0.0141458 -89.37 7.53e-16 ***
## Perc_OC -0.0086994 0.0003263 -26.66 1.27e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0311 on 10 degrees of freedom
## Multiple R-squared: 0.9861, Adjusted R-squared: 0.9847
## F-statistic: 710.6 on 1 and 10 DF, p-value: 1.274e-10 plot(percoc.lm)
ggplot(percoc.40, aes(x=Perc_OC, y=box_prop_released)) +
geom_smooth(method = "lm") +
geom_point(shape=16, size=1) +
theme_bw() ## `geom_smooth()` using formula 'y ~ x'
percoc.lm <- lm(box_prop_released ~ Perc_OC, percoc.40)
summary(percoc.lm)##
## Call:
## lm(formula = box_prop_released ~ Perc_OC, data = percoc.40)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.047503 -0.021611 0.006348 0.018300 0.044623
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.2641718 0.0141458 -89.37 7.53e-16 ***
## Perc_OC -0.0086994 0.0003263 -26.66 1.27e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0311 on 10 degrees of freedom
## Multiple R-squared: 0.9861, Adjusted R-squared: 0.9847
## F-statistic: 710.6 on 1 and 10 DF, p-value: 1.274e-10 plot(percoc.lm)