Hessian fly damage to wheat varieties

Format

block: block factor, 4 levels
genotype factor, 16 wheat varieties
lat: latitude, numeric
long: longitude, numeric
y: number of damaged plants
n: number of total plants

Details

The response is binomial.

Each plot was square.

Source

C. A. Gotway and W. W. Stroup. A Generalized Linear Model Approach to Spatial Data Analysis and Prediction Journal of Agricultural, Biological, and Environmental Statistics, 2, 157-178.

https://doi.org/10.2307/1400401

References

The GLIMMIX procedure. https://www.ats.ucla.edu/stat/SAS/glimmix.pdf

Examples

# \dontrun{

  library(agridat)
  data(gotway.hessianfly)
  dat <- gotway.hessianfly
  
  dat$prop <- dat$y / dat$n
  
  libs(desplot)
  desplot(dat, prop~long*lat,
          aspect=1, # true aspect
          out1=block, num=gen, cex=.75,
          main="gotway.hessianfly")
  

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

  # spaMM package example
  libs(spaMM)
#> Registered S3 methods overwritten by 'registry':
#>   method               from 
#>   print.registry_field proxy
#>   print.registry_entry proxy
#> spaMM (Rousset & Ferdy, 2014, version 3.5.0) is loaded.
#> Type 'help(spaMM)' for a short introduction,
#> 'news(package='spaMM')' for news,
#> and 'citation(spaMM)' for proper citation.
#> 
#> Attaching package: 'spaMM'
#> The following object is masked from 'package:mgcv':
#> 
#>     negbin
  m1 = HLCor(cbind(y, n-y) ~ 1 + gen + (1|block) + Matern(1|long+lat),
             data=dat, family=binomial(), ranPars=list(nu=0.5, rho=1/.7))
  summary(m1)
#> formula: cbind(y, n - y) ~ 1 + gen + (1 | block) + Matern(1 | long + lat)
#> Estimation of lambda by Laplace REML approximation (p_bv).
#> Estimation of fixed effects by Laplace ML approximation (p_v).
#> Family: binomial ( link = logit ) 
#>  ------------ Fixed effects (beta) ------------
#>             Estimate Cond. SE t-value
#> (Intercept)   2.1958   0.6424  3.4181
#> genG02       -0.5166   0.8296 -0.6227
#> genG03       -0.9161   0.8380 -1.0932
#> genG04       -1.8245   0.8010 -2.2777
#> genG05       -0.7466   0.8636 -0.8645
#> genG06       -1.5772   0.8311 -1.8977
#> genG07       -1.1878   0.8337 -1.4248
#> genG08       -2.2601   0.8212 -2.7523
#> genG09       -2.0346   0.7992 -2.5457
#> genG10       -1.2166   0.8583 -1.4175
#> genG11       -2.3692   0.8228 -2.8794
#> genG12       -1.9480   0.8070 -2.4139
#> genG13       -4.6014   0.9629 -4.7785
#> genG14       -3.2362   0.8525 -3.7960
#> genG15       -3.0618   0.8279 -3.6984
#> genG16       -4.1374   0.8887 -4.6554
#>  --------------- Random effects ---------------
#> Family: gaussian ( link = identity ) 
#>                    --- Correlation parameters:
#>    2.rho     2.nu 
#> 1.428571 0.500000 
#>            --- Variance parameters ('lambda'):
#> lambda = var(u) for u ~ Gaussian; 
#>    block  :  4.807e-08 
#>    long + lat  :  0.8619  
#>              --- Coefficients for log(lambda):
#>       Group        Term Estimate Cond.SE
#>       block (Intercept)   -16.85   932.1
#>  long + lat (Intercept)  -0.1486  0.2781
#> # of obs: 64; # of groups: block, 4; long + lat, 64 
#>  ------------- Likelihood values  -------------
#>                        logLik
#> p_v(h) (marginal L): -129.594
#>   p_beta,v(h) (ReL): -123.561
  fixef(m1)
#> (Intercept)      genG02      genG03      genG04      genG05      genG06 
#>   2.1957920  -0.5166378  -0.9160873  -1.8245162  -0.7466366  -1.5771913 
#>      genG07      genG08      genG09      genG10      genG11      genG12 
#>  -1.1877620  -2.2600861  -2.0345601  -1.2165977  -2.3692284  -1.9480056 
#>      genG13      genG14      genG15      genG16 
#>  -4.6014047  -3.2361680  -3.0618020  -4.1374257 
  # The following line fails with "Invalid graphics state"
  # when trying to use pkgdown::build_site
  # filled.mapMM(m1)

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

  # Block random.  See Glimmix manual, output 1.18.
  # Note: (Different parameterization)
  
  libs(lme4)
  l2 <- glmer(cbind(y, n-y) ~ gen + (1|block), data=dat, family=binomial,
    control=glmerControl(check.nlev.gtr.1="ignore"))
  coef(l2)
#> $block
#>    (Intercept)     genG02     genG03    genG04     genG05     genG06     genG07
#> B1    1.507497 -0.1938602 -0.5408074 -1.434189 -0.2037025 -0.9783249 -0.6040858
#> B2    1.502923 -0.1938602 -0.5408074 -1.434189 -0.2037025 -0.9783249 -0.6040858
#> B3    1.493663 -0.1938602 -0.5408074 -1.434189 -0.2037025 -0.9783249 -0.6040858
#> B4    1.509754 -0.1938602 -0.5408074 -1.434189 -0.2037025 -0.9783249 -0.6040858
#>       genG08    genG09     genG10    genG11   genG12   genG13    genG14
#> B1 -1.677434 -1.398435 -0.6817111 -1.462973 -1.45909 -3.55288 -2.507322
#> B2 -1.677434 -1.398435 -0.6817111 -1.462973 -1.45909 -3.55288 -2.507322
#> B3 -1.677434 -1.398435 -0.6817111 -1.462973 -1.45909 -3.55288 -2.507322
#> B4 -1.677434 -1.398435 -0.6817111 -1.462973 -1.45909 -3.55288 -2.507322
#>       genG15    genG16
#> B1 -2.087175 -2.969689
#> B2 -2.087175 -2.969689
#> B3 -2.087175 -2.969689
#> B4 -2.087175 -2.969689
#> 
#> attr(,"class")
#> [1] "coef.mer"

# }