Reading the data

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

Analyzig the data

Definition of model

We will use the following allometric function: \[ W = a L^b \]

Simple linear model

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)

Parameter estimates

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

Confidence intervals

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

Diagnostics

Showing the results

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)