Forage: Forage Effect on Dry Matter Intake

In byandell/pda: Practial Data Analysis for Designed Experiments

Description

The data for this problem are abstracted from a research project at the US Center for Dairy Forage. The experiment involved comparing the effects of five different diets on dry matter intake (dmi), the amount of food eaten. The animals under study were heifers (cows with their first offspring) and mature cows during the lactation (milk-producing) cycle following birth of a calf. The five diets ranged from low (diet=1) to high (diet=5) alfalfa content. The experiment was actually fairly complicated. The proportion of alfalfa for the low group began at 45% and increased to 65% over course of the lactaction cycle, while the high began at 85% and increased to 95%. However, the order of the groups (low to high) remained the same. Interest here focuses on dry matter intake for cows.

Usage

data(ForCow)
data(Forage)

Format

ForCow data frame with 87 observations on 5 variables.

[,1]	hc	factor	Heifer or Cow
[,2]	trt	factor	treatment group
[,3]	code	factor	plot code
[,4]	id	factor	cow identifier
[,5]	dmi	numeric	dry matter intake (DMI)

Forage data frame with 261 observations on 6 variables.

[,1]	cow	factor	cow identifier
[,2]	hc	factor	Heifer or Cow
[,3]	trt	factor	treatment group
[,4]	per	factor	period (Early/Middle/Late)
[,5]	dmi	numeric	dry matter intake (DMI)
[,6]	coce	factor	plot code

Source

Tilak Dhiman

References

Dhiman TR, Kleinmans J, Tessmann NJ, Radloff HD and LD Satter (1995) Digestion and energy balance in lactating dairy cows fed varying ratios of alfalfa silage and grain, J Dairy Science 78, 330.

Dhiman TR and Satter LD (1996) Rumen degradable protein and its effect on microbial protein synthesis, J Dairy Science 00.

Examples

# Average measurement per cow
data( ForCow )

# trt is ordered factor
ForCow$trt <- ordered( ForCow$trt )

# I do heifer, you do cow
ForCowh <- ForCow[ForCow$hc=="H", ]

# stem-and-leaf plot of raw data
stem( ForCowh$dmi, 5 )
# box-plots by treatment
print( bwplot( dmi ~ trt, ForCowh ),
  main = "Forage Raw Data by trt" )

# analysis of variance for heifer
ForCowh.aov <- aov( dmi ~ trt, ForCowh )
ForCowh.aov
summary( ForCowh.aov )

# linear model for heifer -- same idea
ForCowh.lm <- lm( dmi ~ trt, ForCowh )
ForCowh.lm
summary( ForCowh.lm )		# linear models summary
summary.aov( ForCowh.lm )		# anova-type summary

# ordinary means
tapply( ForCowh$dmi, ForCowh$trt, mean )
# predicted values for all responses
data.frame( trt=ForCowh$trt, mean=predict( ForCowh.aov ) )
# predicted values and their standard errors by trt
ForCowh.mean <- predict( ForCowh.aov, se=TRUE )
ForCowh.mean <- data.frame( trt=ForCowh$trt, mean=ForCowh.mean$fit,
   se=ForCowh.mean$se.fit )
ForCowh.mean[ !duplicated( ForCowh$trt ), ]
# least squares means
lsmean( ForCowh.aov )

# polynomial contrasts
coef( ForCowh.aov ) / sqrt(c(1,10,14,10,70))	# or coef( ForCowh.lm )
# use summary.lm to get standard errors
summary.lm( ForCowh.aov )		# or summary( ForCowh.lm )

# S-Plus plot of anova object
# plot( ForCowh.aov )

# residual plot with symbols
tmpdata = data.frame( x = predict( ForCowh.aov ),
  y = resid( ForCowh.aov ), group=ForCowh$trt )
print( xyplot( y ~ x, tmpdata, groups = group,
  xlab="mean by trt", ylab="residual" ))
# residual plot with jittered symbols
# add horizontal lines for 0 and +/- 1 SD
tmpdata$x = jitter( tmpdata$x )
print( xyplot( y ~ x, tmpdata, groups = group,
  panel = function(...) {
    panel.superpose(...)
    panel.abline( h = 0, lty = 2 )
    panel.abline( h=c(-1,1) * std.dev( ForCowh.aov ), lty=3 )
  },
  xlab="mean by trt", ylab="residual" ))

# Full set of cow measurements (3 periods per cow )
data( Forage )

Forage$code <- factor( Forage$code )
Forage$trt <- ordered( Forage$trt )
Forage$per <- ordered( Forage$per, c("Early","Middle","Late") )
Forage$ldmi <- log10( Forage$dmi )

Forage.fit <- aov( ldmi~hc*trt*per+Error( cow ),
   Forage, qr = TRUE )
summary( Forage.fit )

# projections
Forage.proj <- proj( Forage.fit )
Forage.wpint <- apply( Forage.proj$cow[,
   c("hc","trt","hc:trt") ], 1, sum )
Forage.spint <- apply( Forage.proj$Within[,
   c("per","hc:per","trt:per","hc:trt:per") ], 1, sum )
Forage.mean <- Forage.proj[[1]][1]

# H:23.1 Forage Treatment by Period (Heifer or Cow)

# Separate Analysis by Period

ylabs <- c("dry matter intake (DMI)","","")
xlabs <- c("(a) early period","(b) middle period","(c) late period")
mains = c("","Figure H:23.1","")
yaxis <- seq( 12, 24, 2 )
ylim <- log10( c(11.5,24.5) )
lper = levels( Forage$per )

for ( i in seq( along = lper )) {
   tmpfor <- Forage[Forage$per== lper[i], ]
   fit <- lm( log10( dmi )~trt*hc, tmpfor )
   ci.plot( fit, type="b",
     panel = function(x,y,...) {
       panel.abline( h=Forage.mean, lty=3 )
       panel.superpose(x,y,...)
     },
     xlab=xlabs[i], ylab=ylabs[i], main = mains[i],
     more = ( i != 3 ), split=c(i,1,3,1))
} 
#   axis( 2, log10( yaxis ), yaxis )

byandell/pda documentation built on May 13, 2019, 9:27 a.m.