Skip to content

Commit

Permalink
use if(require('asreml')) in examples
Browse files Browse the repository at this point in the history
  • Loading branch information
kwstat committed Feb 2, 2024
1 parent 0f9b997 commit 70615bf
Show file tree
Hide file tree
Showing 72 changed files with 1,635 additions and 1,625 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: agridat
Title: Agricultural Datasets
Version: 1.22
Version: 1.23
Authors@R: person("Kevin","Wright", email="[email protected]",
comment=c(ORCID = "0000-0002-0617-8673"),
role=c("aut","cre","cph"))
Expand Down Expand Up @@ -62,4 +62,4 @@ BugReports: https://github.com/kwstat/agridat/issues
VignetteBuilder: knitr
Language: en-US
Encoding: UTF-8
RoxygenNote: 7.2.3
RoxygenNote: 7.3.1
9 changes: 7 additions & 2 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -1,13 +1,18 @@
# agridat 1.23
# agridat 1.23 (2024-01-30)

## New datasets added

## Other notes

* Removed suggested use of sf/terra packages in examples.

* Add `if(require("asreml", quietly=TRUE))` to example sections. This way we can allow people (or GitHub Actions) to force run the `dontrun` sections, even if `asreml` is not installed.

* cochran.eelworms - Fix a typo reported by U.Genschel, and added more columns for grain yield, straw yield, weeds. Updates to documentation.

* gartner.corn - Remove 'rgdal' package from example (Issue #11).


# agridat 1.22 (2023-08-24)

## New datasets added
Expand Down Expand Up @@ -189,7 +194,7 @@ beaven.barley, perry.springwheat, ridout.appleshoots



# agridat 1.9 - (2014-07-02)
# agridat 1.9 (2014-07-02)

## New datasets added

Expand Down
21 changes: 21 additions & 0 deletions cran-comments.md
Original file line number Diff line number Diff line change
@@ -1,3 +1,24 @@
# agridat 1.23

* Removed suggested use of sf/terra packages in examples.

## Test environments and results

1. R 4.3.2 on Windows 10, devtools::check(cran=TRUE)
2. Rhub Windows Server 2022 R-devel
3. WinBuilder R-Devel

OK (Except usual random notes from Rhub)

## revdepcheck results

* We checked 6 reverse dependencies that were OK.
* The geneticae package failed, but the error message indicates that the failure is not caused by the 'agridat' package.



# -----

# agridat 1.22

* Switched the package to MIT license.
Expand Down
24 changes: 12 additions & 12 deletions man/bailey.cotton.uniformity.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -69,26 +69,26 @@

# Data check. Matches Bailey 1926 Table 1. 28.13, , 46.02, 31.74, 13.52
libs(dplyr)
dat %>% group_by(env) %>% dplyr::summarize(mn=mean(yield))
# dat %>% group_by(env) %>% dplyr::summarize(mn=mean(yield))

libs(desplot)
desplot(dat, yield ~ col*row|env)
desplot(dat, yield ~ col*row|env, main="bailey.cotton.uniformity")

# The yield scales are quite different at each loc, and the dimensions
# are different, so plot each location separately.
# Note: Bailey does not say if plots are 7x15 meters, or 15x7 meters.
# The choices here seem most likely in our opinion.
filter(dat, env=="1921 Sakha") %>%
desplot(yield ~ col*row, main="1921 Sakha", aspect=(20*8.5)/(8*15))
filter(dat, env=="1921 Gemmeiza") %>%
desplot(yield ~ col*row, main="1921 Gemmeiza", aspect=(20*8.5)/(8*15))
filter(dat, env=="1922 Gemmeiza") %>%
desplot(yield ~ col*row, main="1922 Gemmeiza", aspect=(20*8.5)/(8*15))
filter(dat, env=="1921 Giza") %>%
desplot(yield ~ col*row, main="1921 Giza", aspect=(11*6)/(14*8.5))
desplot(dat, yield ~ col*row, subset= env=="1921 Sakha",
main="1921 Sakha", aspect=(20*8.5)/(8*15))
desplot(dat, yield ~ col*row, subset= env=="1921 Gemmeiza",
main="1921 Gemmeiza", aspect=(20*8.5)/(8*15))
desplot(dat, yield ~ col*row, subset= env=="1922 Gemmeiza",
main="1922 Gemmeiza", aspect=(20*8.5)/(8*15))
desplot(dat, yield ~ col*row, subset= env=="1921 Giza",
main="1921 Giza", aspect=(11*6)/(14*8.5))
# 1923 Giza has alternately hi/lo yield rows. Not noticed by Bailey.
filter(dat, env=="1923 Giza") %>%
desplot(yield ~ col*row, main="1923 Giza", aspect=(20*6)/(8*8.5))
desplot(dat, yield ~ col*row, subset= env=="1923 Giza",
main="1923 Giza", aspect=(20*6)/(8*8.5))

}
}
Expand Down
48 changes: 25 additions & 23 deletions man/barrero.maize.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -60,31 +60,33 @@
scales=list(x=list(rot=90)))

# Table 6 of Barrero. Model equation 1.
pacman::p_load(dplyr, asreml, lucid)
dat <- arrange(dat, env)
dat <- mutate(dat,
yearf=factor(year), env=factor(env),
loc=factor(loc), gen=factor(gen), rep=factor(rep))
if(require("asreml", quietly=TRUE)){
libs(dplyr,lucid)
dat <- arrange(dat, env)
dat <- mutate(dat,
yearf=factor(year), env=factor(env),
loc=factor(loc), gen=factor(gen), rep=factor(rep))

m1 <- asreml(yield ~ loc + yearf + loc:yearf, data=dat,
random = ~ gen + rep:loc:yearf +
gen:yearf + gen:loc +
gen:loc:yearf,
residual = ~ dsum( ~ units|env) )

# Variance components for yield match Barrero table 6.
lucid::vc(m1)[1:5,]
## effect component std.error z.ratio bound %ch
## rep:loc:yearf 0.111 0.01092 10 P 0
## gen 0.505 0.03988 13 P 0
## gen:yearf 0.05157 0.01472 3.5 P 0
## gen:loc 0.02283 0.0152 1.5 P 0.2
## gen:loc:yearf 0.2068 0.01806 11 P 0

summary(vc(m1)[6:112,"component"]) # Means match last row of table 6
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1286 0.3577 0.5571 0.8330 1.0322 2.9867
m1 <- asreml(yield ~ loc + yearf + loc:yearf, data=dat,
random = ~ gen + rep:loc:yearf +
gen:yearf + gen:loc +
gen:loc:yearf,
residual = ~ dsum( ~ units|env),
workspace="500mb")

# Variance components for yield match Barrero table 6.
lucid::vc(m1)[1:5,]
## effect component std.error z.ratio bound %ch
## rep:loc:yearf 0.111 0.01092 10 P 0
## gen 0.505 0.03988 13 P 0
## gen:yearf 0.05157 0.01472 3.5 P 0
## gen:loc 0.02283 0.0152 1.5 P 0.2
## gen:loc:yearf 0.2068 0.01806 11 P 0

summary(vc(m1)[6:112,"component"]) # Means match last row of table 6
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1286 0.3577 0.5571 0.8330 1.0322 2.9867
}
}
}
\keyword{datasets}
2 changes: 1 addition & 1 deletion man/battese.survey.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@
Battese, George E and Harter, Rachel M and Fuller, Wayne A. (1988).
An error-components model for prediction of county crop areas using
survey and satellite data.
emph{Journal of the American Statistical Association}, 83, 28-36.
Journal of the American Statistical Association, 83, 28-36.
https://doi.org/10.2307/2288915

Battese (1982) preprint version.
Expand Down
30 changes: 16 additions & 14 deletions man/belamkar.augmented.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -58,21 +58,23 @@
main="belamkar.augmented")

# Belamkar supplement S3 has R code for analysis
library(asreml)
if(require("asreml", quietly=TRUE)){
library(asreml)

# AR1xAR1 model to calculate BLUEs for a single loc
d1 <- droplevels(subset(dat, loc=="Lincoln"))
d1$colf <- factor(d1$col)
d1$rowf <- factor(d1$row)
d1$gen <- factor(d1$gen)
d1$gen_check <- factor(d1$gen_check)
d1 <- d1[order(d1$col),]
d1 <- as.data.frame(d1)
m1 <- asreml(fixed=yield ~ gen_check, data=d1,
random = ~ gen_new:gen,
residual = ~ar1(colf):ar1v(rowf) )
p1 <- predict(m1, classify="gen")
head(p1$pvals)
# AR1xAR1 model to calculate BLUEs for a single loc
d1 <- droplevels(subset(dat, loc=="Lincoln"))
d1$colf <- factor(d1$col)
d1$rowf <- factor(d1$row)
d1$gen <- factor(d1$gen)
d1$gen_check <- factor(d1$gen_check)
d1 <- d1[order(d1$col),]
d1 <- as.data.frame(d1)
m1 <- asreml(fixed=yield ~ gen_check, data=d1,
random = ~ gen_new:gen,
residual = ~ar1(colf):ar1v(rowf) )
p1 <- predict(m1, classify="gen")
head(p1$pvals)
}
}
}
\keyword{datasets}
45 changes: 23 additions & 22 deletions man/besag.bayesian.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -49,30 +49,31 @@
xyplot(yield ~ rrow|col, dat, layout=c(1,3), type='s',
xlab="row", ylab="yield", main="besag.bayesian")

libs(asreml)
if(require("asreml", quietly=TRUE)) {
libs(asreml, lucid)

# Use asreml to fit a model with AR1 gradient in rows
dat <- transform(dat, cf=factor(col), rf=factor(rrow))
m1 <- asreml(yield ~ -1 + gen, data=dat, random= ~ rf)
m1 <- update(m1, random= ~ ar1v(rf))
m1 <- update(m1)
m1 <- update(m1)
m1 <- update(m1)
lucid::vc(m1)
# Use asreml to fit a model with AR1 gradient in rows
dat <- transform(dat, cf=factor(col), rf=factor(rrow))
m1 <- asreml(yield ~ -1 + gen, data=dat, random= ~ rf)
m1 <- update(m1, random= ~ ar1v(rf))
m1 <- update(m1)
m1 <- update(m1)
m1 <- update(m1)
lucid::vc(m1)

# Visualize trends, similar to Besag figure 2.
# Need 'as.vector' because asreml uses a named vector
dat$res <- unname(m1$resid)
dat$geneff <- coef(m1)$fixed[as.numeric(dat$gen)]
dat <- transform(dat, fert=yield-geneff-res)
libs(lattice)
xyplot(geneff ~ rrow|col, dat, layout=c(1,3), type='s',
main="besag.bayesian - Variety effects", ylim=c(5,15 ))
xyplot(fert ~ rrow|col, dat, layout=c(1,3), type='s',
main="besag.bayesian - Fertility", ylim=c(-2,2))
xyplot(res ~ rrow|col, dat, layout=c(1,3), type='s',
main="besag.bayesian - Residuals", ylim=c(-4,4))

# Visualize trends, similar to Besag figure 2.
# Need 'as.vector' because asreml uses a named vector
dat$res <- unname(m1$resid)
dat$geneff <- coef(m1)$fixed[as.numeric(dat$gen)]
dat <- transform(dat, fert=yield-geneff-res)
libs(lattice)
xyplot(geneff ~ rrow|col, dat, layout=c(1,3), type='s',
main="besag.bayesian - Variety effects", ylim=c(5,15 ))
xyplot(fert ~ rrow|col, dat, layout=c(1,3), type='s',
main="besag.bayesian - Fertility", ylim=c(-2,2))
xyplot(res ~ rrow|col, dat, layout=c(1,3), type='s',
main="besag.bayesian - Residuals", ylim=c(-4,4))
}
}
}

Expand Down
40 changes: 20 additions & 20 deletions man/besag.endive.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -74,27 +74,27 @@
## Xy -0.144800 NA NA NA
## eta 0.806200 NA NA NA


libs(asreml)

# Now try an AR1xAR1 model.
dat2 <- transform(dat, xf=factor(col), yf=factor(row),
pres=as.numeric(disease=="Y"))

m2 <- asreml(pres ~ 1, data=dat2,
resid = ~ar1(xf):ar1(yf))
# The 0/1 response is arbitrary, but there is some suggestion
# of auto-correlation in the x (.17) and y (.10) directions,
# suggesting the pattern is more 'patchy' than just random noise,
# but is it meaningful?

libs(lucid)
vc(m2)
## effect component std.error z.ratio bound %ch
## xf:yf(R) 0.1301 0.003798 34 P 0
## xf:yf!xf!cor 0.1699 0.01942 8.7 U 0
## xf:yf!yf!cor 0.09842 0.02038 4.8 U 0

if(require("asreml", quietly=TRUE)) {
libs(asreml,lucid)

# Now try an AR1xAR1 model.
dat2 <- transform(dat, xf=factor(col), yf=factor(row),
pres=as.numeric(disease=="Y"))

m2 <- asreml(pres ~ 1, data=dat2,
resid = ~ar1(xf):ar1(yf))
# The 0/1 response is arbitrary, but there is some suggestion
# of auto-correlation in the x (.17) and y (.10) directions,
# suggesting the pattern is more 'patchy' than just random noise,
# but is it meaningful?

lucid::vc(m2)
## effect component std.error z.ratio bound %ch
## xf:yf(R) 0.1301 0.003798 34 P 0
## xf:yf!xf!cor 0.1699 0.01942 8.7 U 0
## xf:yf!yf!cor 0.09842 0.02038 4.8 U 0
}

}
}
Expand Down
Loading

0 comments on commit 70615bf

Please sign in to comment.