crowder.seeds: Germination of Orobanche seeds for two genotypes and two...

Description Format Details Source References Examples

Description

Number of Orobanche seeds tested/germinated for two genotypes and two treatments.

Format

plate

Factor for replication

gen

Factor for genotype with levels O73, O75

extract

Factor for extract from bean, cucumber

germ

Number of seeds that germinated

n

Total number of seeds tested

Details

Egyptian broomrape, orobanche aegyptiaca is a parasitic plant family. The plants have no chlorophyll and grow on the roots of other plants. The seeds remain dormant in soil until certain compounds from living plants stimulate germination.

Two genotypes were studied in the experiment, O. aegyptiaca 73 and O. aegyptiaca 75. The seeds were brushed with one of two extracts prepared from either a bean plant or cucmber plant.

The experimental design was a 2x2 factorial, each with 5 or 6 reps of plates.

Source

Crowder, M.J., 1978. Beta-binomial anova for proportions. Appl. Statist., 27, 34-37. http://doi.org/10.2307/2346223

References

N. E. Breslow and D. G. Clayton. 1993. Approximate inference in generalized linear mixed models. Journal of the American Statistical Association, 88:9-25. http://doi.org/10.2307/2290687

Y. Lee and J. A. Nelder. 1996. Hierarchical generalized linear models with discussion. J. R. Statist. Soc. B, 58:619-678.

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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
library(agridat)
data(crowder.seeds)
dat <- crowder.seeds
m1.glm <- m1.glmm <- m1.bb <- m1.hglm <- NA


# ----- Graphic
require(lattice)
dotplot(germ/n~gen|extract, dat, main="crowder.seeds")


# ----- GLM.
# family=binomial() fixes dispersion at 1
# family=quasibinomial() estimates dispersion, had larger std errors
m1.glm <- glm(cbind(germ,n-germ) ~ gen*extract,
              data=dat,
              #family="binomial",
              family=quasibinomial()
              )
summary(m1.glm)


# --- GLMM.  Assumes Gaussian random effects
require(MASS)
m1.glmm <- glmmPQL(cbind(germ, n-germ) ~ gen*extract, random= ~1|plate,
  family=binomial(), data=dat)
summary(m1.glmm)


# ----- AODS3 package
# require(aods3)
# m1.bb <- aodml(cbind(germ, n-germ) ~ gen * extract, data=dat, family="bb")


# ----- HGML package. Beta-binomial with beta-distributed random effects
# require(hglm)
# m1.hglm <- hglm(fixed= germ/n ~ I(gen=="O75")*extract, weights=n, data=dat,
#                 random=~1|plate, family=binomial(), rand.family=Beta(),
#                 fix.disp=1)

## Not run: 
# Compare coefficients

round(summary(m1.glm)$coef,2)
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)               -0.41       0.25   -1.64     0.12
## genO75                    -0.15       0.30   -0.48     0.64
## extractcucumber            0.54       0.34    1.58     0.13
## genO75:extractcucumber     0.78       0.42    1.86     0.08

round(summary(m1.glmm)$tTable,2)
##                        Value Std.Error DF t-value p-value
## (Intercept)            -0.44      0.25 17   -1.80    0.09
## genO75                 -0.10      0.31 17   -0.34    0.74
## extractcucumber         0.52      0.34 17    1.56    0.14
## genO75:extractcucumber  0.80      0.42 17    1.88    0.08

## round(summary(m1.bb)$BCoef,2)
##                        Estimate Std. Error z value Pr(> |z|)
## (Intercept)               -0.44       0.22   -2.04      0.04
## genO75                    -0.10       0.27   -0.36      0.72
## extractcucumber            0.52       0.30    1.76      0.08
## genO75:extractcucumber     0.80       0.38    2.11      0.03

## round(summary(m1.hglm)$FixCoefMat,2)
##                                     Estimate Std. Error t-value Pr(>|t|)
## (Intercept)                            -0.47       0.24   -1.92     0.08
## I(gen == "O75")TRUE                    -0.08       0.31   -0.25     0.81
## extractcucumber                         0.51       0.33    1.53     0.16
## I(gen == "O75")TRUE:extractcucumber     0.83       0.43    1.92     0.08

## End(Not run)

## Not run: 

  # --- rjags version ---

# JAGS/BUGS.  See http://mathstat.helsinki.fi/openbugs/Examples/Seeds.html
# Germination rate depends on p, which is a logit of a linear predictor
# based on genotype and extract, plus random deviation to intercept

# To match the output on the BUGS web page, use: dat$gen=="O73".
# We use dat$gen=="O75" to compare with the parameterization above.
jdat =list(germ = dat$germ, n = dat$n,
           root = as.numeric(dat$extract=="cucumber"),
           gen = as.numeric(dat$gen=="O75"),
           nobs = nrow(dat))

jinit = list(int = 0, genO75 = 0, extcuke = 0, g75ecuke = 0, tau = 10)

# Use logical names (unlike BUGS documentation)
mod.bug =
"model {
  for(i in 1:nobs) {
    germ[i] ~ dbin(p[i], n[i])
    b[i] ~ dnorm(0.0, tau)
    logit(p[i]) <- int + genO75 * gen[i] + extcuke * root[i] +
                   g75ecuke * gen[i] * root[i] + b[i]
  }
  int ~ dnorm(0.0, 1.0E-6)
  genO75 ~ dnorm(0.0, 1.0E-6)
  extcuke ~ dnorm(0.0, 1.0E-6)
  g75ecuke ~ dnorm(0.0, 1.0E-6)
  tau ~ dgamma(0.001, 0.001)
  sigma <- 1 / sqrt(tau)
}"

require(rjags)
oo <- textConnection(mod.bug)
j1 <- jags.model(oo, data=jdat, inits=jinit, n.chains=1)
close(oo)

c1 <- coda.samples(j1, c("int","genO75","g75ecuke","extcuke","sigma"),
                   n.iter=20000)
summary(c1) # Medians are very similar to estimates from hglm
# require(lucid)
# print(vc(c1),3)
##             Mean    SD    2.5
## extcuke   0.543  0.331 -0.118   0.542   1.2   
## g75ecuke  0.807  0.436 -0.0586  0.802   1.7   
## genO75   -0.0715 0.309 -0.665  -0.0806  0.581 
## int      -0.479  0.241 -0.984  -0.473  -0.0299
## sigma     0.289  0.142  0.0505  0.279   0.596 

# Plot observed data with HPD intervals for germination probability
c2 <- coda.samples(j1, c("p"), n.iter=20000)
hpd <- HPDinterval(c2)[[1]]
med <- summary(c2, quantiles=.5)$quantiles
fit <- data.frame(med, hpd)

if(require(latticeExtra)){
  obs <- dotplot(1:21 ~ germ/n, dat,
                 main="crowder.seeds", ylab="plate",
                 col=as.numeric(dat$gen), pch=substring(dat$extract,1))
  obs + segplot(1:21 ~ lower + upper, data=fit, centers=med)
}

# ----------------------------------------------------------------------------
## --- R2jags version ---

require("agridat")
require("R2jags")
dat <- crowder.seeds

# To match the output on the BUGS web page, use: dat$gen=="O73".
# We use dat$gen=="O75" to compare with the parameterization above.
jdat =list(germ = dat$germ, n = dat$n,
           root = as.numeric(dat$extract=="cucumber"),
           gen = as.numeric(dat$gen=="O75"),
           nobs = nrow(dat))

jinit = list(list(int = 0, genO75 = 0, extcuke = 0, g75ecuke = 0, tau = 10))

mod.bug = function() {
  for(i in 1:nobs) {
    germ[i] ~ dbin(p[i], n[i])
    b[i] ~ dnorm(0.0, tau)
    logit(p[i]) <- int + genO75 * gen[i] + extcuke * root[i] +
      g75ecuke * gen[i] * root[i] + b[i]
  }
  int ~ dnorm(0.0, 1.0E-6)
  genO75 ~ dnorm(0.0, 1.0E-6)
  extcuke ~ dnorm(0.0, 1.0E-6)
  g75ecuke ~ dnorm(0.0, 1.0E-6)
  tau ~ dgamma(0.001, 0.001)
  sigma <- 1 / sqrt(tau)
}

parms <- c("int","genO75","g75ecuke","extcuke","sigma")

j1 <- jags(data=jdat, inits=jinit, parms, model.file=mod.bug,
           n.iter=20000, n.chains=1)
print(j1)
##          mu.vect sd.vect   2.5
## extcuke    0.519   0.325 -0.140  0.325   0.531   0.728   1.158
## g75ecuke   0.834   0.429 -0.019  0.552   0.821   1.101   1.710
## genO75    -0.096   0.305 -0.670 -0.295  -0.115   0.089   0.552
## int       -0.461   0.236 -0.965 -0.603  -0.455  -0.312   0.016
## sigma      0.255   0.148  0.033  0.140   0.240   0.352   0.572
## deviance 103.319   7.489 90.019 98.010 102.770 108.689 117.288

traceplot(as.mcmc(j1))
densityplot(as.mcmc(j1))
HPDinterval(as.mcmc(j1))


## End(Not run)

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