Turfgrass ratings for different treatments
Turfgrass ratings for different treatments
A data frame with 128 observations on the following 6 variables.
weekweek number
repblocking factor
managemanagement factor, 4 levels
nitronitrogen factor, 2 levels
ratingturfgrass rating, 4 ordered levels
countnumber of samples for a given rating
Turf color was assessed on a scale of Poor, Average, Good, Excellent.
The data are the number of times that a combination of management style and nitrogen level received a particular rating across four replicates and four sampling weeks. The eight treatments were in a completely randomized design.
Nitrogen level 1 is 2.5 g/m^2, level 2 is 5 g/m^2.
Management 1 = N applied with no supplemental water injection.
M2 = surface applied with supplemental water injection.
M3 = nitrogen injected 7.6 cm deep
M4 = nitrogen injected 12.7 cm deep.
Schabenberger, Oliver and Francis J. Pierce. 2002. Contemporary Statistical Models for the Plant and Soil Sciences. CRC Press. Page 380.
library(agridat)
data(karcher.turfgrass)
dat <- karcher.turfgrass
dat$rating <- ordered(dat$rating, levels=c('Poor','Average', 'Good','Excellent'))
ftable(xtabs(~manage+nitro+rating, dat)) # Table 6.19 of Schabenberger
# Probably would choose management M3, nitro N2
mosaicplot(xtabs(count ~ manage + rating + nitro, dat),
shade=TRUE, dir=c('h','v','h'),
main="karcher.turfgrass - turfgrass ratings")
# Multinomial logistic model. Probit Ordered Logistic Regression.
libs(MASS)
m1 <- polr(rating ~ nitro*manage + week, dat, weights=count, Hess=TRUE, method='logistic')
summary(m1)
# Try to match the "predicted marginal probability distribution" of
# Schabenberger table 6.20. He doesn't define "marginal".
# Are the interaction terms included before aggregation?
# Are 'margins' calculated before/after back-transforming?
# At what level is the covariate 'week' included?
# Here is what Schabenberger presents:
## M1 M2 M3 M4 | N1 N2
## Poor .668 .827 .001 .004 | .279 .020
## Avg .330 .172 .297 .525 | .712 .826
## Good .002 .001 .695 .008 | .008 .153
## Exc .000 .000 .007 .003 | .001 .001
## We use week=3.5, include interactions, then average
newd <- expand.grid(manage=levels(dat$manage), nitro=levels(dat$nitro), week=3.5)
newd <- cbind(newd, predict(m1, newdata=newd, type='probs')) # probs)
print(aggregate( . ~ manage, data=newd, mean), digits=2)
## manage nitro week Poor Average Good Excellent
## 1 M1 1.5 3.5 0.67 0.33 0.0011 0.0000023
## 2 M2 1.5 3.5 0.76 0.24 0.00059 0.0000012
## 3 M3 1.5 3.5 0.0023 0.48 0.52 0.0042
## 4 M4 1.5 3.5 0.0086 0.57 0.42 0.0035Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.