Number of lambs born to 3 breeds on 3 farms

data("mead.lamb")

Format

A data frame with 36 observations on the following 4 variables.

farm: farm: F1, F2, F3
breed: breed: B1, B2, B3
lambclass: lambing class: L0, L1, L2, L3
y: count of ewes in class

Details

The data 'y' are counts of ewes in different lambing classes. The classes are number of live lambs per birth for 0, 1, 2, 3+ lambs.

Source

Roger Mead, Robert N Curnow, Anne M Hasted. 2002. Statistical Methods in Agriculture and Experimental Biology, 3rd ed. Chapman and Hall. Page 359.

References

None

Examples

# \dontrun{

library(agridat)
data(mead.lamb)
dat <- mead.lamb

# farm 1 has more ewes in lambclass 3
d2 <- xtabs(y ~ farm+breed+lambclass, data=dat)
mosaicplot(d2, color=c("lemonchiffon1","moccasin","lightsalmon1","indianred"),
           xlab="farm/lambclass", ylab="breed", main="mead.lamb")

names(dat) <- c('F','B','L','y') # for compactness
# Match totals in Mead example 14.6
libs(dplyr)
dat <- group_by(dat, F,B)
summarize(dat, y=sum(y))
#> `summarise()` regrouping output by 'F' (override with `.groups` argument)
#> # A tibble: 9 x 3
#> # Groups:   F [3]
#>   F     B         y
#>   <fct> <fct> <int>
#> 1 F1    A       150
#> 2 F1    B        46
#> 3 F1    C        78
#> 4 F2    A        72
#> 5 F2    B        79
#> 6 F2    C        28
#> 7 F3    A       224
#> 8 F3    B       129
#> 9 F3    C        34
##          F      B     y
##   <fctr> <fctr> <int>
## 1     F1      A   150
## 2     F1      B    46
## 3     F1      C    78
## 4     F2      A    72
## 5     F2      B    79
## 6     F2      C    28
## 7     F3      A   224
## 8     F3      B   129
## 9     F3      C    34


# Models
m1 <- glm(y ~ F + B + F:B, data=dat,
          family=poisson(link=log))
m2 <- update(m1, y ~ F + B + F:B + L)
m3 <- update(m1, y ~ F + B + F:B + L + B:L)
m4 <- update(m1, y ~ F + B + F:B + L + F:L)
m5 <- update(m1, y ~ F + B + F:B + L + B:L + F:L)

AIC(m1, m2, m3, m4, m5) # Model 4 has best AIC
#>    df      AIC
#> m1  9 852.9800
#> m2 12 306.5457
#> m3 18 303.5781
#> m4 18 206.1520
#> m5 24 213.8873
##    df      AIC
## m1  9 852.9800
## m2 12 306.5457
## m3 18 303.5781
## m4 18 206.1520
## m5 24 213.8873

# Change contrasts for Miroslav
m4 <- update(m4,
             contrasts=list(F=contr.sum,B=contr.sum,L=contr.sum))
summary(m4)
#> 
#> Call:
#> glm(formula = y ~ F + B + L + F:B + F:L, family = poisson(link = log), 
#>     data = dat, contrasts = list(F = contr.sum, B = contr.sum, 
#>         L = contr.sum))
#> 
#> Deviance Residuals: 
#>      Min        1Q    Median        3Q       Max  
#> -2.00000  -0.49709   0.04627   0.47793   1.54534  
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  2.38679    0.07147  33.395  < 2e-16 ***
#> F1           0.23967    0.08795   2.725 0.006431 ** 
#> F2          -0.44980    0.11371  -3.956 7.63e-05 ***
#> B1           0.56941    0.05214  10.922  < 2e-16 ***
#> B2           0.02240    0.05790   0.387 0.698885    
#> L1          -0.50321    0.11242  -4.476 7.60e-06 ***
#> L2           0.45373    0.08930   5.081 3.76e-07 ***
#> L3           1.42046    0.07553  18.806  < 2e-16 ***
#> F1:B1        0.04256    0.07061   0.603 0.546678    
#> F2:B1       -0.28552    0.08104  -3.523 0.000426 ***
#> F1:B2       -0.59242    0.08541  -6.936 4.03e-12 ***
#> F2:B2        0.35428    0.08405   4.215 2.50e-05 ***
#> F1:L1       -0.34198    0.15261  -2.241 0.025031 *  
#> F2:L1        0.01066    0.17867   0.060 0.952420    
#> F1:L2       -0.76829    0.12326  -6.233 4.58e-10 ***
#> F2:L2        0.12824    0.13935   0.920 0.357432    
#> F1:L3       -0.05738    0.09454  -0.607 0.543896    
#> F2:L3        0.23542    0.12030   1.957 0.050351 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 1008.825  on 35  degrees of freedom
#> Residual deviance:   18.845  on 18  degrees of freedom
#> AIC: 206.15
#> 
#> Number of Fisher Scoring iterations: 5
#> 

# Match deviance table from Mead
libs(broom)
all <- do.call(rbind, lapply(list(m1, m2, m3, m4, m5), broom::glance))
all$model <- unlist(lapply(list(m1, m2, m3, m4, m5),
                               function(x) as.character(formula(x)[3])))
all[,c('model','deviance','df.residual')]
#> # A tibble: 5 x 3
#>   model                       deviance df.residual
#>   <chr>                          <dbl>       <int>
#> 1 F + B + F:B                    684.           27
#> 2 F + B + L + F:B                131.           24
#> 3 F + B + L + F:B + B:L          116.           18
#> 4 F + B + L + F:B + F:L           18.8          18
#> 5 F + B + L + F:B + B:L + F:L     14.6          12
##                         model  deviance df.residual
## 1                 F + B + F:B 683.67257          27
## 2             F + B + L + F:B 131.23828          24
## 3       F + B + L + F:B + B:L 116.27069          18
## 4       F + B + L + F:B + F:L  18.84460          18
## 5 F + B + L + F:B + B:L + F:L  14.57987          12


if(0){
  # Using MASS::loglm
  libs(MASS)
  # Note: without 'fitted=TRUE', devtools::run_examples has an error
  m4b <- MASS::loglm(y ~ F + B + F:B + L + F:L, data = dat, fitted=TRUE)
  # Table of farm * class interactions. Match Mead p. 360
  round(coef(m4b)$F.L,2)
  fitted(m4b)
  resid(m4b)
  # libs(vcd)
  # mosaic(m4b, shade=TRUE,
  # formula = ~ F + B + F:B + L + F:L,
  # residual_type="rstandard", keep_aspect=FALSE)
}

# }