kenward.cattle: Repeated measurement of weights of calves with two...

Description Usage Format Details Source References Examples

Description

Repeated measurements of the weights of calves from a trial on the control of intestinal parasites.

Usage

1
data("kenward.cattle")

Format

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

animal

animal factor

trt

treatment factor, A or B

day

day, numberic, 0-133

weight

bodyweight, kg

Details

Grazing cattle can ingest larvae, which deprives the host animal of nutrients and weakens the immune system, affecting the growth of the animal.

Two treatments A and B were applied randomly to 60 animals (30 each in two groups) to control the disease.

Each animal was weighed 11 times at two-week intervals (one week between the final two measurements).

Is there a difference in treatments, and when does that difference first become manifest?

Source

Kenward, Michael G. (1987). A Method for Comparing Profiles of Repeated Measurements. Applied Statistics, 36, 296-308. Table 1. http://doi.org/10.2307/2347788

References

W. Zhang, C. Leng and C. Y. Tang (2015). A joint modelling approach for longitudinal studies J. R. Statist. Soc. B, 77 (2015), 219–238. http://doi.org/10.1111/rssb.12065

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
data(kenward.cattle)
dat <- kenward.cattle

# Profile plots
require(lattice)
foo1 <- xyplot(weight~day|trt, data=dat, type='l', group=animal,
               xlab="Day", ylab="Animal weight", main="kenward.cattle")
print(foo1)

# ----------------------------------------------------------------------------

## Not run: 

  # lme4. Fixed treatment intercepts, treatment polynomial trend.
  # Random deviation for each animal
  require(lme4)
  m1a <-lmer(weight ~ trt*poly(day, 4) + (1|animal), data=dat,
             REML = FALSE)
  # Change separate polynomials into common polynomial
  m1b <-lmer(weight ~ trt + poly(day, 4) + (1|animal), data=dat,
             REML = FALSE)
  # Drop treatment differences
  m1c <-lmer(weight ~ poly(day, 4) + (1|animal), data=dat,
             REML = FALSE)
  anova(m1a, m1b, m1c) # Significant differences between trt polynomials

  # Overlay polynomial predictions on plot
  require(latticeExtra)
  dat$pred <- predict(m1a, re.form=NA)
  foo1 + xyplot(pred ~ day|trt, data=dat,
                lwd=2, col="black", type='l')

  # A Kenward-Roger Approximation and Parametric Bootstrap
  library(pbkrtest)
  KRmodcomp(m1b, m1c) # Non-signif
  # Model comparison of nested models using parametric bootstrap methods
  # PBmodcomp(m1b, m1c, nsim=500)
  ## Parametric bootstrap test; time: 13.20 sec; samples: 500 extremes: 326;
  ## large : weight ~ trt + poly(day, 4) + (1 | animal)
  ## small : weight ~ poly(day, 4) + (1 | animal)
  ##          stat df p.value
  ## LRT    0.2047  1  0.6509
  ## PBtest 0.2047     0.6527


## End(Not run)

# ----------------------------------------------------------------------------

## Not run: 
  # ASREML approach to model. Not final by any means.
  # Maybe a spline curve for each treatment, plus random deviations for each time
  # asreml3
  require(asreml)
  m1 <- asreml(weight ~  1 + lin(day) +    # overall line
                 trt + trt:lin(day),       # different line for each treatment
               data=dat,
               random = ~ spl(day) +       # overall spline
                 trt:spl(day) +            # different spline for each treatment
                 dev(day) + trt:dev(day) ) # non-spline deviation at each time*trt
  
  p1 <- predict(m1, data=dat, classify="trt:day")
  p1 <- p1$predictions$pvals
  foo2 <- xyplot(predicted.value ~ day|trt, p1, type='l', lwd=2, lty=1, col="black")
  
  require(latticeExtra)
  print(foo1 + foo2)

  # Not much evidence for treatment differences
  
  anova(m1)
  ##               Df Sum of Sq Wald statistic Pr(Chisq)    
  ## (Intercept)    1  37128459         139060    <2e-16 ***
  ## trt            1       455              2    0.1917    
  ## lin(day)       1    570798           2138    <2e-16 ***
  ## trt:lin(day)   1       283              1    0.3031    
  ## residual (MS)          267                             
  
  require(lucid)
  vc(m1)
  ##               effect component std.error z.ratio constr
  ##             spl(day)  25.29    24.09        1       pos
  ##             dev(day)   1.902    4.923       0.39    pos
  ## trt:spl(day)!trt.var   0.00003  0.000002   18      bnd 
  ## trt:dev(day)!trt.var   0.00003  0.000002   18      bnd 
  ##           R!variance 267       14.84       18       pos

## End(Not run)

# ----------------------------------------------------------------------------

## Not run: 
  # ASREML approach to model. Not final by any means.
  # Maybe a spline curve for each treatment, plus random deviations for each time
  ## require(asreml4)
  ## m1 <- asreml(weight ~  1 + lin(day) +    # overall line
  ##                trt + trt:lin(day),       # different line for each treatment
  ##              data=dat,
  ##              random = ~ spl(day) +       # overall spline
  ##                trt:spl(day) +            # different spline for each treatment
  ##                dev(day) + trt:dev(day) ) # non-spline deviation at each time*trt
  
  ## p1 <- predict(m1, data=dat, classify="trt:day",
  ##               design.points=list(day=0:133))
  ## p1 <- p1$pvals
  ## foo2 <- xyplot(predicted.value ~ day|trt, p1, type='l', lwd=2, lty=1, col="black")
  
  ## require(latticeExtra)
  ## print(foo1 + foo2)

  ## # Not much evidence for treatment differences
  
  ## wald(m1)
  ## ##               Df Sum of Sq Wald statistic Pr(Chisq)    
  ## ## (Intercept)    1  37128459         139060    <2e-16 ***
  ## ## trt            1       455              2    0.1917    
  ## ## lin(day)       1    570798           2138    <2e-16 ***
  ## ## trt:lin(day)   1       283              1    0.3031    
  ## ## residual (MS)          267                             
  
  ## require(lucid)
  ## vc(m1)
  ## ##               effect component std.error z.ratio constr
  ## ##             spl(day)  25.29    24.09        1       pos
  ## ##             dev(day)   1.902    4.923       0.39    pos
  ## ## trt:spl(day)!trt.var   0.00003  0.000002   18      bnd 
  ## ## trt:dev(day)!trt.var   0.00003  0.000002   18      bnd 
  ## ##           R!variance 267       14.84       18       pos


## End(Not run)

# ----------------------------------------------------------------------------

agridat documentation built on May 2, 2019, 4:01 p.m.