Data <- read.csv("Rainbowtrout.csv")
names(Data)
## [1] "Length" "Weight"
head(Data, 10)
## Length Weight
## 1 7.5 10
## 2 9.2 9
## 3 9.5 25
## 4 9.7 35
## 5 10.5 20
## 6 10.6 35
## 7 19.0 70
## 8 19.0 100
## 9 19.1 70
## 10 19.5 68
tail(Data, 5)
## Length Weight
## 35 44.5 970
## 36 45.6 980
## 37 48.7 1700
## 38 49.2 1290
## 39 51.7 1650
Length <- Data$Length
Weight <- Data$Weight
We will use the following allometric function: \[ W = a L^b \]
We then conduct a simple linear regression analysis with the following simple linear model: \[ \log W = \log a + b \log L. \]
res <- lm(log(Weight)~log(Length))
summary(res)
##
## Call:
## lm(formula = log(Weight) ~ log(Length))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59889 -0.08977 -0.00832 0.14493 0.69146
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.04898 0.28442 -10.72 6.67e-13 ***
## log(Length) 2.60235 0.08474 30.71 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.286 on 37 degrees of freedom
## Multiple R-squared: 0.9622, Adjusted R-squared: 0.9612
## F-statistic: 943.1 on 1 and 37 DF, p-value: < 2.2e-16
plot(log(Length), log(Weight), xlim=c(0,4), ylim=c(0,8))
abline(res)
names(res)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
coef(res)
## (Intercept) log(Length)
## -3.048982 2.602353
alpha.est <- as.numeric(coef(res)[1])
beta.est <- as.numeric(coef(res)[2])
a.est <- exp(alpha.est)
b.est <- beta.est
Est <- c(a.est, b.est)
Est
## [1] 0.04740718 2.60235273
confint(res)
## 2.5 % 97.5 %
## (Intercept) -3.625272 -2.472691
## log(Length) 2.430651 2.774054
alpha.ci <- confint(res)[1,]
beta.ci <- confint(res)[2,]
#alpha.ci
#beta.ci
a.ci <- exp(alpha.ci)
b.ci <- beta.ci
CI <- rbind(a.ci, b.ci)
#CI
Out <- cbind(Est, CI)
row.names(Out) <- c("a", "b")
Out
## Est 2.5 % 97.5 %
## a 0.04740718 0.02664184 0.08435757
## b 2.60235273 2.43065137 2.77405408
plot(Length, Weight, xlab="Length (mm)",
ylab="Weight (g)", main="With CI for estimates")
pred <- exp(predict(res, int="c", level = 0.95))
matlines(Length, pred, lty=c(1,2,2), col=1)
plot(Length, Weight, xlab="Length (mm)",
ylab="Weight (g)", main="With CI for estimates")
pred <- exp(predict(res, int="p", level = 0.95))
## Warning in predict.lm(res, int = "p", level = 0.95): predictions on current data refer to _future_ responses
matlines(Length, pred, lty=c(1,2,2), col=1)