## Initialization

> library(NCStats)

# Salmon Sperm Example

## Background

Vladic et al. (2002) recorded (in SalmonSperm.csv) the probability of successful egg fertilization (fert.success) and the length of sperm tail end piece (step.len). They asked “Are fertilization success and length of sperm related?”

> setwd("C:/aaaWork/Web/GitHub/NCMTH207/modules/SLRegression")
> ss <- read.csv("SalmonSperm.csv")
> ss <- ss[-c(1,10,11),]  # only for class demo purposes
> str(ss)
'data.frame':   11 obs. of  3 variables:
$step.len : num 2.94 3 3.02 3.17 3.18 3.2 3.27 3.31 3.72 3.84 ...$ fert.succ: num  3 2.2 7 7 13.5 10.4 6.7 12.8 37.8 30 ...
$mat : Factor w/ 2 levels "Adult","Parr": 2 2 1 2 1 1 1 1 2 2 ... > xlbl <- "Sperm Tail End Piece Length (um)" > ylbl <- "Fertilization Success" ## Lecture Support I – Model Fitting and Simple Predictions > ( lm1 <- lm(fert.succ~step.len,data=ss) ) Coefficients: (Intercept) step.len -100.21 34.61  > fitPlot(lm1,xlab=xlbl,ylab=ylbl) > predict(lm1,data.frame(step.len=3.5))  1 20.92912  ## Lecture Support II – Sampling Variability > summary(lm1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -100.205 13.015 -7.699 3.00e-05 step.len 34.610 3.889 8.901 9.35e-06 Residual standard error: 4.366 on 9 degrees of freedom Multiple R-squared: 0.898, Adjusted R-squared: 0.8866 F-statistic: 79.22 on 1 and 9 DF, p-value: 9.35e-06  > confint(lm1)  2.5 % 97.5 % (Intercept) -129.64815 -70.76202 step.len 25.81336 43.40619 > fitPlot(lm1,interval="both",xlab=xlbl,ylab=ylbl) > predict(lm1,data.frame(step.len=3.5),interval="confidence")  fit lwr upr 1 20.92912 17.5967 24.26153 > predict(lm1,data.frame(step.len=3.5),interval="prediction")  fit lwr upr 1 20.92912 10.50502 31.35321 > predictionPlot(lm1,data.frame(step.len=c(3.3,3.5)),interval="prediction", xlab=xlbl,ylab=ylbl)  obs step.len fit lwr upr 1 1 3.3 14.00716 3.687506 24.32682 2 2 3.5 20.92912 10.505016 31.35321 ## Lecture Support III – Model Comparisons > anova(lm1) Analysis of Variance Table Response: fert.succ Df Sum Sq Mean Sq F value Pr(>F) step.len 1 1510.23 1510.23 79.219 9.35e-06 Residuals 9 171.58 19.06  ## Lecture Support IV – Assumption Checking > residPlot(lm1) > adTest(lm1$residuals)
Anderson-Darling normality test with x
A = 0.4022, p-value = 0.2962
> outlierTest(lm1)

No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferonni p
12 3.717896          0.0058892     0.064781

# Petrels Example

## Background

Croxall (1982) examined the weight loss of adult petrels during periods of egg incubation. He examined 13 species but some had measurements for both sexes such that 19 measurements are found in Petrels.csv. For each measurement the mean initial weight (g) and mean weight lost (g/g/d) were recorded. Determine if the mean initial weight significant explains variability in mean weight lost.

> petrels <- read.csv("Petrels.csv")
> str(petrels)
'data.frame':   19 obs. of  4 variables:
$species : Factor w/ 13 levels "Diomedea chrysostoma",..: 2 2 4 4 1 1 3 3 3 9 ...$ sex        : Factor w/ 4 levels "both","female",..: 3 2 3 2 3 2 3 2 1 3 ...
$weight : int 10577 9022 3922 3694 3751 3624 3305 3000 2996 668 ...$ weight.loss: num  0.0087 0.0096 0.013 0.011 0.01 0.012 0.011 0.0125 0.0109 0.0128 ...

## Analysis

> lm1 <- lm(weight.loss~weight,data=petrels)
> fitPlot(lm1,xlab="Weight (g)",ylab="Weight Loss (g/g/d)")

> residPlot(lm1)

> with(petrels,max(weight)/min(weight))
[1] 251.8333
> ## transChooser(lm1) # interactive, results not shown
> petrels$log.wt <- log(petrels$weight)
> petrels$log.wtloss <- log(petrels$weight.loss)
> lm2 <- lm(log.wtloss~log.wt,data=petrels)
> fitPlot(lm2,xlab="log Weight (g)",ylab="log Weight Loss (g/g/d)")

> residPlot(lm2)

> adTest(lm2\$residuals)
Anderson-Darling normality test with x
A = 0.3881, p-value = 0.3514
> anova(lm2)
Analysis of Variance Table

Response: log.wtloss
Df Sum Sq Mean Sq F value    Pr(>F)
log.wt     1 6.5113  6.5113  140.65 1.204e-09
Residuals 17 0.7870  0.0463                  
> summary(lm2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.73403    0.19792  -8.761 1.04e-07
log.wt      -0.33632    0.02836 -11.860 1.20e-09

Residual standard error: 0.2152 on 17 degrees of freedom
Multiple R-squared: 0.8922, Adjusted R-squared: 0.8858
F-statistic: 140.6 on 1 and 17 DF,  p-value: 1.204e-09 
> confint(lm2)
                 2.5 %     97.5 %
(Intercept) -2.1516113 -1.3164546
log.wt      -0.3961507 -0.2764885
> ( p.log.wtloss <- predict(lm2,data.frame(log.wt=log(5000)),interval="confidence") )
        fit       lwr       upr
1 -4.598532 -4.746569 -4.450495
> exp(p.log.wtloss)*exp(anova(lm2)[2,3]/2)
         fit         lwr        upr
1 0.01030234 0.008884726 0.01194614