From 70615bf91b553a93d0df2b5076d21a672c7657df Mon Sep 17 00:00:00 2001 From: Kevin Wright Date: Fri, 2 Feb 2024 11:26:04 -0600 Subject: [PATCH] use if(require('asreml')) in examples --- DESCRIPTION | 4 +- NEWS.md | 9 +- cran-comments.md | 21 +++ man/bailey.cotton.uniformity.Rd | 24 ++-- man/barrero.maize.Rd | 48 +++---- man/battese.survey.Rd | 2 +- man/belamkar.augmented.Rd | 30 +++-- man/besag.bayesian.Rd | 45 +++---- man/besag.endive.Rd | 40 +++--- man/besag.met.Rd | 154 +++++++++++----------- man/besag.triticale.Rd | 52 ++++---- man/box.cork.Rd | 66 +++++----- man/bridges.cucumber.Rd | 81 ++++++------ man/burgueno.alpha.Rd | 132 +++++++++---------- man/burgueno.rowcol.Rd | 43 +++--- man/burgueno.unreplicated.Rd | 62 ++++----- man/butron.maize.Rd | 14 +- man/cleveland.soil.Rd | 52 ++++---- man/cochran.eelworms.Rd | 48 +++---- man/cullis.earlygen.Rd | 91 ++++++------- man/damesa.maize.Rd | 11 +- man/digby.jointregression.Rd | 8 +- man/diggle.cow.Rd | 32 ++--- man/draper.safflower.uniformity.Rd | 9 +- man/durban.rowcol.Rd | 48 +++---- man/durban.splitplot.Rd | 41 +++--- man/federer.diagcheck.Rd | 78 +++++------ man/fisher.barley.Rd | 75 +++++------ man/gartner.corn.Rd | 146 ++++++++++----------- man/gilmour.serpentine.Rd | 85 ++++++------ man/gilmour.slatehall.Rd | 44 ++++--- man/gomez.multilocsplitplot.Rd | 4 +- man/hadasch.lettuce.Rd | 32 ++--- man/hanks.sprinkler.Rd | 69 +++------- man/hansen.multi.uniformity.Rd | 29 +++-- man/harris.wateruse.Rd | 188 +++++++++++++-------------- man/hayman.tobacco.Rd | 6 +- man/heady.fertilizer.Rd | 2 +- man/hildebrand.systems.Rd | 76 +++++------ man/ilri.sheep.Rd | 40 +++--- man/john.alpha.Rd | 63 +++++---- man/kalamkar.wheat.uniformity.Rd | 26 ++-- man/kempton.barley.uniformity.Rd | 24 ++-- man/kempton.slatehall.Rd | 35 +++-- man/kenward.cattle.Rd | 60 ++++----- man/lonnquist.maize.Rd | 35 ++--- man/payne.wheat.Rd | 87 +++++++------ man/perry.springwheat.Rd | 30 +++-- man/piepho.barley.uniformity.Rd | 31 ++--- man/piepho.cocksfoot.Rd | 47 +++---- man/shafi.tomato.uniformity.Rd | 2 +- man/sinclair.clover.Rd | 2 +- man/smith.corn.uniformity.Rd | 2 +- man/smith.wheat.uniformity.Rd | 1 + man/snedecor.asparagus.Rd | 110 ++++++++-------- man/steptoe.morex.pheno.Rd | 48 +++---- man/stroup.nin.Rd | 54 ++++---- man/stroup.splitplot.Rd | 70 +++++----- man/tesfaye.millet.Rd | 49 +++---- man/thompson.cornsoy.Rd | 4 +- man/usgs.herbicides.Rd | 3 +- man/verbyla.lupin.Rd | 154 +++++++++++----------- man/vold.longterm.Rd | 3 +- man/vsn.lupin3.Rd | 168 ++++++++++++------------ man/wedderburn.barley.Rd | 5 +- man/weiss.incblock.Rd | 25 ++-- man/weiss.lattice.Rd | 27 ++-- man/wiebe.wheat.uniformity.Rd | 1 + man/yates.oats.Rd | 47 ++----- vignettes/agridat_data.Rmd | 2 +- vignettes/agridat_graphical_gems.Rmd | 2 +- vignettes/agridat_intro.Rmd | 2 +- 72 files changed, 1635 insertions(+), 1625 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5c946bd6..984ca8dc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: agridat Title: Agricultural Datasets -Version: 1.22 +Version: 1.23 Authors@R: person("Kevin","Wright", email="kw.stat@gmail.com", comment=c(ORCID = "0000-0002-0617-8673"), role=c("aut","cre","cph")) @@ -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 diff --git a/NEWS.md b/NEWS.md index 833c21eb..9762e27c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -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 @@ -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 diff --git a/cran-comments.md b/cran-comments.md index 03aca2f2..fd790a2b 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -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. diff --git a/man/bailey.cotton.uniformity.Rd b/man/bailey.cotton.uniformity.Rd index 2d138843..aa3ff620 100644 --- a/man/bailey.cotton.uniformity.Rd +++ b/man/bailey.cotton.uniformity.Rd @@ -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)) } } diff --git a/man/barrero.maize.Rd b/man/barrero.maize.Rd index ad1c0d49..12e4c5ba 100644 --- a/man/barrero.maize.Rd +++ b/man/barrero.maize.Rd @@ -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} diff --git a/man/battese.survey.Rd b/man/battese.survey.Rd index 52e8a020..adb932f6 100644 --- a/man/battese.survey.Rd +++ b/man/battese.survey.Rd @@ -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. diff --git a/man/belamkar.augmented.Rd b/man/belamkar.augmented.Rd index 3a75e7b2..a3d6ec44 100644 --- a/man/belamkar.augmented.Rd +++ b/man/belamkar.augmented.Rd @@ -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} diff --git a/man/besag.bayesian.Rd b/man/besag.bayesian.Rd index e0469a3e..a34f9c51 100644 --- a/man/besag.bayesian.Rd +++ b/man/besag.bayesian.Rd @@ -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)) + } } } diff --git a/man/besag.endive.Rd b/man/besag.endive.Rd index f0224fa5..02090e76 100644 --- a/man/besag.endive.Rd +++ b/man/besag.endive.Rd @@ -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 + } } } diff --git a/man/besag.met.Rd b/man/besag.met.Rd index 2198e874..5e8f72d3 100644 --- a/man/besag.met.Rd +++ b/man/besag.met.Rd @@ -85,85 +85,85 @@ ## 91.90 210.75 63.03 112.05 28.39 237.36 72.72 42.97 ## ... etc ... - - libs(asreml) - - # Average reps - datm <- aggregate(yield ~ county + gen, data=dat, FUN=mean) - # asreml Using 'rcov' ALWAYS requires sorting the data - datm <- datm[order(datm$gen),] - - m1 <- asreml(yield ~ gen, data=datm, - random = ~ county, - residual = ~ dsum( ~ units|gen)) - libs(lucid) - vc(m1)[1:7,] - ## effect component std.error z.ratio bound %ch - ## county 1324 836.1 1.6 P 0.2 - ## gen_G01!R 91.98 58.91 1.6 P 0.1 - ## gen_G02!R 210.6 133.6 1.6 P 0.1 - ## gen_G03!R 63.06 40.58 1.6 P 0.1 - ## gen_G04!R 112.1 71.59 1.6 P 0.1 - ## gen_G05!R 28.35 18.57 1.5 P 0.2 - ## gen_G06!R 237.4 150.8 1.6 P 0 + if(require("asreml", quietly=TRUE)) { + libs(asreml, lucid) + + # Average reps + datm <- aggregate(yield ~ county + gen, data=dat, FUN=mean) + # asreml Using 'rcov' ALWAYS requires sorting the data + datm <- datm[order(datm$gen),] + + m1 <- asreml(yield ~ gen, data=datm, + random = ~ county, + residual = ~ dsum( ~ units|gen)) + vc(m1)[1:7,] + ## effect component std.error z.ratio bound %ch + ## county 1324 836.1 1.6 P 0.2 + ## gen_G01!R 91.98 58.91 1.6 P 0.1 + ## gen_G02!R 210.6 133.6 1.6 P 0.1 + ## gen_G03!R 63.06 40.58 1.6 P 0.1 + ## gen_G04!R 112.1 71.59 1.6 P 0.1 + ## gen_G05!R 28.35 18.57 1.5 P 0.2 + ## gen_G06!R 237.4 150.8 1.6 P 0 - # We get the same results from asreml & lme - # plot(m1$vparameters[-1], - # m1l$sigma^2 * c(1, coef(m1l$modelStruct$varStruct, unc = FALSE))^2) + # We get the same results from asreml & lme + # plot(m1$vparameters[-1], + # m1l$sigma^2 * c(1, coef(m1l$modelStruct$varStruct, unc = FALSE))^2) + + # The following example shows how to construct a GxE biplot + # from the FA2 model. + + + dat <- besag.met + dat <- transform(dat, xf=factor(col), yf=factor(row)) + dat <- dat[order(dat$county, dat$xf, dat$yf), ] + + # First, AR1xAR1 + m1 <- asreml(yield ~ county, data=dat, + random = ~ gen:county, + residual = ~ dsum( ~ ar1(xf):ar1(yf)|county)) + # Add FA1 + m2 <- update(m1, random=~gen:fa(county,1)) # rotate.FA=FALSE + # FA2 + m3 <- update(m2, random=~gen:fa(county,2)) + asreml.options(extra=50) + m3 <- update(m3, maxit=50) + asreml.options(extra=0) + + # Use the loadings to make a biplot + vars <- vc(m3) + psi <- vars[grepl("!var$", vars$effect), "component"] + la1 <- vars[grepl("!fa1$", vars$effect), "component"] + la2 <- vars[grepl("!fa2$", vars$effect), "component"] + mat <- as.matrix(data.frame(psi, la1, la2)) + # I tried using rotate.fa=FALSE, but it did not seem to + # give orthogonal vectors. Rotate by hand. + rot <- svd(mat[,-1])$v # rotation matrix + lam <- mat[,-1] %*% rot # Rotate the loadings + colnames(lam) <- c("load1", "load2") + + co3 <- coef(m3)$random # Scores are the GxE coefficients + ix1 <- grepl("_Comp1$", rownames(co3)) + ix2 <- grepl("_Comp2$", rownames(co3)) + sco <- matrix(c(co3[ix1], co3[ix2]), ncol=2, byrow=FALSE) + sco <- sco %*% rot # Rotate the scores + dimnames(sco) <- list(levels(dat$gen) , c('load1','load2')) + rownames(lam) <- levels(dat$county) + sco[,1:2] <- -1 * sco[,1:2] + lam[,1:2] <- -1 * lam[,1:2] + biplot(sco, lam, cex=.5, main="FA2 coefficient biplot (asreml)") + # G variance matrix + gvar <- lam %*% t(lam) + diag(mat[,1]) - # The following example shows how to construct a GxE biplot - # from the FA2 model. - - - dat <- besag.met - dat <- transform(dat, xf=factor(col), yf=factor(row)) - dat <- dat[order(dat$county, dat$xf, dat$yf), ] - - # First, AR1xAR1 - m1 <- asreml(yield ~ county, data=dat, - random = ~ gen:county, - residual = ~ dsum( ~ ar1(xf):ar1(yf)|county)) - # Add FA1 - m2 <- update(m1, random=~gen:fa(county,1)) # rotate.FA=FALSE - # FA2 - m3 <- update(m2, random=~gen:fa(county,2)) - asreml.options(extra=50) - m3 <- update(m3, maxit=50) - asreml.options(extra=0) - - # Use the loadings to make a biplot - vars <- vc(m3) - psi <- vars[grepl("!var$", vars$effect), "component"] - la1 <- vars[grepl("!fa1$", vars$effect), "component"] - la2 <- vars[grepl("!fa2$", vars$effect), "component"] - mat <- as.matrix(data.frame(psi, la1, la2)) - # I tried using rotate.fa=FALSE, but it did not seem to - # give orthogonal vectors. Rotate by hand. - rot <- svd(mat[,-1])$v # rotation matrix - lam <- mat[,-1] %*% rot # Rotate the loadings - colnames(lam) <- c("load1", "load2") - - co3 <- coef(m3)$random # Scores are the GxE coefficients - ix1 <- grepl("_Comp1$", rownames(co3)) - ix2 <- grepl("_Comp2$", rownames(co3)) - sco <- matrix(c(co3[ix1], co3[ix2]), ncol=2, byrow=FALSE) - sco <- sco %*% rot # Rotate the scores - dimnames(sco) <- list(levels(dat$gen) , c('load1','load2')) - rownames(lam) <- levels(dat$county) - sco[,1:2] <- -1 * sco[,1:2] - lam[,1:2] <- -1 * lam[,1:2] - biplot(sco, lam, cex=.5, main="FA2 coefficient biplot (asreml)") - # G variance matrix - gvar <- lam %*% t(lam) + diag(mat[,1]) - - # Now get predictions and make an ordinary biplot - p3 <- predict(m3, data=dat, classify="county:gen") - p3 <- p3$pvals - libs("gge") - bi3 <- gge(p3, predicted.value ~ gen*county, scale=FALSE) - if(interactive()) dev.new() - # Very similar to the coefficient biplot - biplot(bi3, stand=FALSE, main="SVD biplot of FA2 predictions") + # Now get predictions and make an ordinary biplot + p3 <- predict(m3, data=dat, classify="county:gen") + p3 <- p3$pvals + libs("gge") + bi3 <- gge(p3, predicted.value ~ gen*county, scale=FALSE) + if(interactive()) dev.new() + # Very similar to the coefficient biplot + biplot(bi3, stand=FALSE, main="SVD biplot of FA2 predictions") + } } } diff --git a/man/besag.triticale.Rd b/man/besag.triticale.Rd index 9a9dab2b..6b0cb0dc 100644 --- a/man/besag.triticale.Rd +++ b/man/besag.triticale.Rd @@ -66,32 +66,34 @@ as.table=TRUE, type="s", layout=c(1,3), main="besag.triticale") - libs(asreml) + if(require("asreml", quietly=TRUE)) { + libs(asreml) - # Besag uses an adjustment based on neighboring plots. - # This analysis fits the standard AR1xAR1 residual model - - dat <- dat[order(dat$xf, dat$yf), ] - m3 <- asreml(yield ~ gen + rate + nitro + regulator + - gen:rate + gen:nitro + gen:regulator + - rate:nitro + rate:regulator + - nitro:regulator + yf, data=dat, - resid = ~ ar1(xf):ar1(yf)) - wald(m3) # Strongly significant gen, rate, regulator - ## Df Sum of Sq Wald statistic Pr(Chisq) - ## (Intercept) 1 1288255 103.971 < 2.2e-16 *** - ## gen 2 903262 72.899 < 2.2e-16 *** - ## rate 1 104774 8.456 0.003638 ** - ## nitro 1 282 0.023 0.880139 - ## regulator 2 231403 18.676 8.802e-05 *** - ## yf 2 3788 0.306 0.858263 - ## gen:rate 2 1364 0.110 0.946461 - ## gen:nitro 2 30822 2.488 0.288289 - ## gen:regulator 4 37269 3.008 0.556507 - ## rate:nitro 1 1488 0.120 0.728954 - ## rate:regulator 2 49296 3.979 0.136795 - ## nitro:regulator 2 41019 3.311 0.191042 - ## residual (MS) 12391 + # Besag uses an adjustment based on neighboring plots. + # This analysis fits the standard AR1xAR1 residual model + + dat <- dat[order(dat$xf, dat$yf), ] + m3 <- asreml(yield ~ gen + rate + nitro + regulator + + gen:rate + gen:nitro + gen:regulator + + rate:nitro + rate:regulator + + nitro:regulator + yf, data=dat, + resid = ~ ar1(xf):ar1(yf)) + wald(m3) # Strongly significant gen, rate, regulator + ## Df Sum of Sq Wald statistic Pr(Chisq) + ## (Intercept) 1 1288255 103.971 < 2.2e-16 *** + ## gen 2 903262 72.899 < 2.2e-16 *** + ## rate 1 104774 8.456 0.003638 ** + ## nitro 1 282 0.023 0.880139 + ## regulator 2 231403 18.676 8.802e-05 *** + ## yf 2 3788 0.306 0.858263 + ## gen:rate 2 1364 0.110 0.946461 + ## gen:nitro 2 30822 2.488 0.288289 + ## gen:regulator 4 37269 3.008 0.556507 + ## rate:nitro 1 1488 0.120 0.728954 + ## rate:regulator 2 49296 3.979 0.136795 + ## nitro:regulator 2 41019 3.311 0.191042 + ## residual (MS) 12391 + } } } diff --git a/man/box.cork.Rd b/man/box.cork.Rd index 90621962..b57dbfa4 100644 --- a/man/box.cork.Rd +++ b/man/box.cork.Rd @@ -18,12 +18,10 @@ } \source{ - C.R. Rao (1948). Tests of significance in multivariate analysis. \emph{Biometrika}, 35, 58-79. - https://doi.org/10.2307/2332629 - + https://doi.org/10.2307/2332629 } \references{ @@ -66,42 +64,42 @@ line.col=rep(c("royalblue","red","#009900","dark orange", "#999999","#a6761d","deep pink"), length=nrow(dat2))) - - # asreml - libs(asreml) - # Unstructured covariance - dat$dir <- factor(dat$dir) - dat$tree <- factor(dat$tree) - dat <- dat[order(dat$tree, dat$dir), ] - - # Unstructured covariance matrix - m1 <- asreml(y~dir, data=dat, residual = ~ tree:us(dir)) + if(require("asreml", quietly=TRUE)) { + libs(asreml, lucid) - libs(lucid) - vc(m1) + # Unstructured covariance + dat$dir <- factor(dat$dir) + dat$tree <- factor(dat$tree) + dat <- dat[order(dat$tree, dat$dir), ] - # Note: 'rcor' is a personal function to extract the correlations - # into a matrix format - # round(kw::rcor(m1)$dir, 2) - # E N S W - # E 219.93 223.75 229.06 171.37 - # N 223.75 290.41 288.44 226.27 - # S 229.06 288.44 350.00 259.54 - # W 171.37 226.27 259.54 226.00 + # Unstructured covariance matrix + m1 <- asreml(y~dir, data=dat, residual = ~ tree:us(dir)) - # Note: Wolfinger used a common diagonal variance - - # Factor Analytic with different specific variances - # fixme: does not work with asreml4 - # m2 <- update(m1, residual = ~tree:facv(dir,1)) - # round(kw::rcor(m2)$dir, 2) - # E N S W - # E 219.94 209.46 232.85 182.27 - # N 209.46 290.41 291.82 228.43 - # S 232.85 291.82 349.99 253.94 - # W 182.27 228.43 253.94 225.99 + lucid::vc(m1) + # Note: 'rcor' is a personal function to extract the correlations + # into a matrix format + # round(kw::rcor(m1)$dir, 2) + # E N S W + # E 219.93 223.75 229.06 171.37 + # N 223.75 290.41 288.44 226.27 + # S 229.06 288.44 350.00 259.54 + # W 171.37 226.27 259.54 226.00 + + # Note: Wolfinger used a common diagonal variance + + # Factor Analytic with different specific variances + # fixme: does not work with asreml4 + # m2 <- update(m1, residual = ~tree:facv(dir,1)) + # round(kw::rcor(m2)$dir, 2) + # E N S W + # E 219.94 209.46 232.85 182.27 + # N 209.46 290.41 291.82 228.43 + # S 232.85 291.82 349.99 253.94 + # W 182.27 228.43 253.94 225.99 + } + } } \keyword{datasets} diff --git a/man/bridges.cucumber.Rd b/man/bridges.cucumber.Rd index f7f6861f..41e0dd03 100644 --- a/man/bridges.cucumber.Rd +++ b/man/bridges.cucumber.Rd @@ -66,49 +66,50 @@ bwplot(yield ~ loc|factor(.sample), dat20, main="bridges.cucumber - graphical inference") - libs(asreml) # asreml + if(require("asreml", quietly=TRUE)) { + libs(asreml,lucid) - ## Random row/col/resid. Same as Bridges 1989, p. 147 - m1 <- asreml(yield ~ 1 + gen + loc + loc:gen, - random = ~ rowf:loc + colf:loc, data=dat) + ## Random row/col/resid. Same as Bridges 1989, p. 147 + m1 <- asreml(yield ~ 1 + gen + loc + loc:gen, + random = ~ rowf:loc + colf:loc, data=dat) - libs(lucid) - lucid::vc(m1) - ## effect component std.error z.ratio bound %ch - ## rowf:loc 31.62 23.02 1.4 P 0 - ## colf:loc 18.08 15.32 1.2 P 0 - ## units(R) 31.48 12.85 2.4 P 0 - - ## Random row/col/resid at each loc. Matches p. 147 - m2 <- asreml(yield ~ 1 + gen + loc + loc:gen, - random = ~ at(loc):rowf + at(loc):colf, data=dat, - resid = ~ dsum( ~ units|loc)) - lucid::vc(m2) - ## effect component std.error z.ratio bound %ch - ## at(loc, Clemson):rowf 32.32 36.58 0.88 P 0 - ## at(loc, Tifton):rowf 30.92 28.63 1.1 P 0 - ## at(loc, Clemson):colf 22.55 28.78 0.78 P 0 - ## at(loc, Tifton):colf 13.62 14.59 0.93 P 0 - ## loc_Clemson(R) 46.85 27.05 1.7 P 0 - ## loc_Tifton(R) 16.11 9.299 1.7 P 0 - - predict(m2, data=dat, classify='loc:gen')$pvals - ## loc gen predicted.value std.error status - ## 1 Clemson Dasher 45.6 5.04 Estimable - ## 2 Clemson Guardian 31.6 5.04 Estimable - ## 3 Clemson Poinsett 21.4 5.04 Estimable - ## 4 Clemson Sprint 26 5.04 Estimable - ## 5 Tifton Dasher 50.5 3.89 Estimable - ## 6 Tifton Guardian 38.7 3.89 Estimable - ## 7 Tifton Poinsett 33 3.89 Estimable - ## 8 Tifton Sprint 39.2 3.89 Estimable + lucid::vc(m1) + ## effect component std.error z.ratio bound %ch + ## rowf:loc 31.62 23.02 1.4 P 0 + ## colf:loc 18.08 15.32 1.2 P 0 + ## units(R) 31.48 12.85 2.4 P 0 + + ## Random row/col/resid at each loc. Matches p. 147 + m2 <- asreml(yield ~ 1 + gen + loc + loc:gen, + random = ~ at(loc):rowf + at(loc):colf, data=dat, + resid = ~ dsum( ~ units|loc)) + lucid::vc(m2) + ## effect component std.error z.ratio bound %ch + ## at(loc, Clemson):rowf 32.32 36.58 0.88 P 0 + ## at(loc, Tifton):rowf 30.92 28.63 1.1 P 0 + ## at(loc, Clemson):colf 22.55 28.78 0.78 P 0 + ## at(loc, Tifton):colf 13.62 14.59 0.93 P 0 + ## loc_Clemson(R) 46.85 27.05 1.7 P 0 + ## loc_Tifton(R) 16.11 9.299 1.7 P 0 + + predict(m2, data=dat, classify='loc:gen')$pvals + ## loc gen predicted.value std.error status + ## 1 Clemson Dasher 45.6 5.04 Estimable + ## 2 Clemson Guardian 31.6 5.04 Estimable + ## 3 Clemson Poinsett 21.4 5.04 Estimable + ## 4 Clemson Sprint 26 5.04 Estimable + ## 5 Tifton Dasher 50.5 3.89 Estimable + ## 6 Tifton Guardian 38.7 3.89 Estimable + ## 7 Tifton Poinsett 33 3.89 Estimable + ## 8 Tifton Sprint 39.2 3.89 Estimable + + # Is a heterogeneous model justified? Maybe not. + # m1$loglik + ## -67.35585 + # m2$loglik + ## -66.35621 + } - # Is a heterogeneous model justified? Maybe not. - # m1$loglik - ## -67.35585 - # m2$loglik - ## -66.35621 - } } diff --git a/man/burgueno.alpha.Rd b/man/burgueno.alpha.Rd index 7ec6f7ea..73ab2510 100644 --- a/man/burgueno.alpha.Rd +++ b/man/burgueno.alpha.Rd @@ -62,75 +62,75 @@ ## Residual 133200 365 - libs(asreml) + if(require("asreml", quietly=TRUE)) { + libs(asreml) - dat <- transform(dat, xf=factor(col), yf=factor(row)) - dat <- dat[order(dat$xf, dat$yf),] - - # Sequence of models on page 36 - - m1 <- asreml(yield ~ gen, data=dat) - m1$loglik # -232.13 - - m2 <- asreml(yield ~ gen, data=dat, - random = ~ rep) - m2$loglik # -223.48 - - # Inc Block model - m3 <- asreml(yield ~ gen, data=dat, - random = ~ rep/block) - m3$loglik # -221.42 - m3$coef$fixed # Matches solution on p. 27 - - # AR1xAR1 model - m4 <- asreml(yield ~ 1 + gen, data=dat, - resid = ~ar1(xf):ar1(yf)) - m4$loglik # -221.47 - plot(varioGram(m4), main="burgueno.alpha") # Figure 1 - - m5 <- asreml(yield ~ 1 + gen, data=dat, - random= ~ yf, resid = ~ar1(xf):ar1(yf)) - m5$loglik # -220.07 - - m6 <- asreml(yield ~ 1 + gen + pol(yf,-2), data=dat, - resid = ~ar1(xf):ar1(yf)) - m6$loglik # -204.64 - - m7 <- asreml(yield ~ 1 + gen + lin(yf), data=dat, - random= ~ spl(yf), resid = ~ar1(xf):ar1(yf)) - m7$loglik # -212.51 - - m8 <- asreml(yield ~ 1 + gen + lin(yf), data=dat, - random= ~ spl(yf)) - m8$loglik # -213.91 - - # Polynomial model with predictions - m9 <- asreml(yield ~ 1 + gen + pol(yf,-2) + pol(xf,-2), data=dat, - random= ~ spl(yf), - resid = ~ar1(xf):ar1(yf)) - m9 <- update(m9) - m9$loglik # -191.44 vs -189.61 - p9 <- predict(m9, classify="gen:xf:yf", levels=list(xf=1,yf=1)) - p9 - - m10 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat, - resid = ~ar1(xf):ar1(yf)) - m10$loglik # -211.56 - - m11 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat, - random= ~ spl(yf), - resid = ~ar1(xf):ar1(yf)) - m11$loglik # -208.90 + dat <- transform(dat, xf=factor(col), yf=factor(row)) + dat <- dat[order(dat$xf, dat$yf),] + + # Sequence of models on page 36 of Burgueno + + m1 <- asreml(yield ~ gen, data=dat) + m1$loglik # -232.13 + + m2 <- asreml(yield ~ gen, data=dat, + random = ~ rep) + m2$loglik # -223.48 + + # Inc Block model + m3 <- asreml(yield ~ gen, data=dat, + random = ~ rep/block) + m3$loglik # -221.42 + m3$coef$fixed # Matches solution on p. 27 + + # AR1xAR1 model + m4 <- asreml(yield ~ 1 + gen, data=dat, + resid = ~ar1(xf):ar1(yf)) + m4$loglik # -221.47 + plot(varioGram(m4), main="burgueno.alpha") # Figure 1 + + m5 <- asreml(yield ~ 1 + gen, data=dat, + random= ~ yf, resid = ~ar1(xf):ar1(yf)) + m5$loglik # -220.07 + + m6 <- asreml(yield ~ 1 + gen + pol(yf,-2), data=dat, + resid = ~ar1(xf):ar1(yf)) + m6$loglik # -204.64 + + m7 <- asreml(yield ~ 1 + gen + lin(yf), data=dat, + random= ~ spl(yf), resid = ~ar1(xf):ar1(yf)) + m7$loglik # -212.51 + + m8 <- asreml(yield ~ 1 + gen + lin(yf), data=dat, + random= ~ spl(yf)) + m8$loglik # -213.91 + + # Polynomial model with predictions + m9 <- asreml(yield ~ 1 + gen + pol(yf,-2) + pol(xf,-2), data=dat, + random= ~ spl(yf), + resid = ~ar1(xf):ar1(yf)) + m9 <- update(m9) + m9$loglik # -191.44 vs -189.61 - m12 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat, - random= ~ spl(yf)+spl(xf), - resid = ~ar1(xf):ar1(yf)) - m12$loglik # -206.82 + m10 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat, + resid = ~ar1(xf):ar1(yf)) + m10$loglik # -211.56 + + m11 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat, + random= ~ spl(yf), + resid = ~ar1(xf):ar1(yf)) + m11$loglik # -208.90 + + m12 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat, + random= ~ spl(yf)+spl(xf), + resid = ~ar1(xf):ar1(yf)) + m12$loglik # -206.82 + + m13 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat, + random= ~ spl(yf)+spl(xf)) + m13$loglik # -207.52 + } - m13 <- asreml(yield ~ 1 + gen + lin(yf)+lin(xf), data=dat, - random= ~ spl(yf)+spl(xf)) - m13$loglik # -207.52 - } } diff --git a/man/burgueno.rowcol.Rd b/man/burgueno.rowcol.Rd index 34648c7c..2ec0bb6b 100644 --- a/man/burgueno.rowcol.Rd +++ b/man/burgueno.rowcol.Rd @@ -61,27 +61,28 @@ ## rep (Intercept) 0.1916 0.4378 ## Residual 0.1796 0.4238 - libs(asreml) - - # AR1 x AR1 with linear row/col effects, random spline row/col - dat <- transform(dat, xf=factor(col), yf=factor(row)) - dat <- dat[order(dat$xf,dat$yf),] - m2 <- asreml(yield ~ gen + lin(yf) + lin(xf), data=dat, - random = ~ spl(yf) + spl(xf), - resid = ~ ar1(xf):ar1(yf)) - m2 <- update(m2) # More iterations - - # Scaling of spl components has changed in asreml from old versions - libs(lucid) - vc(m2) # Match Burgueno p. 42 - ## effect component std.error z.ratio bound %ch - ## spl(yf) 0.09077 0.08252 1.1 P 0 - ## spl(xf) 0.08107 0.08209 0.99 P 0 - ## xf:yf(R) 0.1482 0.03119 4.8 P 0 - ## xf:yf!xf!cor 0.1152 0.2269 0.51 U 0.1 - ## xf:yf!yf!cor 0.009467 0.2414 0.039 U 0.9 + if(require("asreml", quietly=TRUE)) { + libs(asreml,lucid) + + # AR1 x AR1 with linear row/col effects, random spline row/col + dat <- transform(dat, xf=factor(col), yf=factor(row)) + dat <- dat[order(dat$xf,dat$yf),] + m2 <- asreml(yield ~ gen + lin(yf) + lin(xf), data=dat, + random = ~ spl(yf) + spl(xf), + resid = ~ ar1(xf):ar1(yf)) + m2 <- update(m2) # More iterations + + # Scaling of spl components has changed in asreml from old versions + lucid::vc(m2) # Match Burgueno p. 42 + ## effect component std.error z.ratio bound %ch + ## spl(yf) 0.09077 0.08252 1.1 P 0 + ## spl(xf) 0.08107 0.08209 0.99 P 0 + ## xf:yf(R) 0.1482 0.03119 4.8 P 0 + ## xf:yf!xf!cor 0.1152 0.2269 0.51 U 0.1 + ## xf:yf!yf!cor 0.009467 0.2414 0.039 U 0.9 + + plot(varioGram(m2), main="burgueno.rowcol") + } - plot(varioGram(m2), main="burgueno.rowcol") - } } diff --git a/man/burgueno.unreplicated.Rd b/man/burgueno.unreplicated.Rd index 172ac8be..7f2d6a75 100644 --- a/man/burgueno.unreplicated.Rd +++ b/man/burgueno.unreplicated.Rd @@ -54,38 +54,40 @@ col=check, num=gen, #text=gen, cex=.3, # aspect unknown main="burgueno.unreplicated") - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)) { + libs(asreml,lucid) - # AR1 x AR1 with random genotypes - dat <- transform(dat, xf=factor(col), yf=factor(row)) - dat <- dat[order(dat$xf,dat$yf),] - m2 <- asreml(yield ~ 1, data=dat, random = ~ gen, - resid = ~ ar1(xf):ar1(yf)) - vc(m2) - ## effect component std.error z.ratio bound %ch - ## gen 0.9122 0.127 7.2 P 0 - ## xf:yf(R) 0.4993 0.05601 8.9 P 0 - ## xf:yf!xf!cor -0.2431 0.09156 -2.7 U 0 - ## xf:yf!yf!cor 0.1255 0.07057 1.8 U 0.1 + # AR1 x AR1 with random genotypes + dat <- transform(dat, xf=factor(col), yf=factor(row)) + dat <- dat[order(dat$xf,dat$yf),] + m2 <- asreml(yield ~ 1, data=dat, random = ~ gen, + resid = ~ ar1(xf):ar1(yf)) + lucid::vc(m2) + ## effect component std.error z.ratio bound %ch + ## gen 0.9122 0.127 7.2 P 0 + ## xf:yf(R) 0.4993 0.05601 8.9 P 0 + ## xf:yf!xf!cor -0.2431 0.09156 -2.7 U 0 + ## xf:yf!yf!cor 0.1255 0.07057 1.8 U 0.1 + + # Note the strong saw-tooth pattern in the variogram. Seems to + # be column effects. + plot(varioGram(m2), xlim=c(0,15), ylim=c(0,9), zlim=c(0,0.5), + main="burgueno.unreplicated - AR1xAR1") + # libs(lattice) # Show how odd columns are high + # bwplot(resid(m2) ~ col, data=dat, horizontal=FALSE) + + # Define an even/odd column factor as fixed effect + # dat$oddcol <- factor(dat$col %% 2) + # The modulus operator throws a bug, so do it the hard way. + dat$oddcol <- factor(dat$col - floor(dat$col / 2) *2 ) - # Note the strong saw-tooth pattern in the variogram. Seems to - # be column effects. - plot(varioGram(m2), xlim=c(0,15), ylim=c(0,9), zlim=c(0,0.5), - main="burgueno.unreplicated - AR1xAR1") - # libs(lattice) # Show how odd columns are high - # bwplot(resid(m2) ~ col, data=dat, horizontal=FALSE) - - # Define an even/odd column factor as fixed effect - # dat$oddcol <- factor(dat$col %% 2) - # The modulus operator throws a bug, so do it the hard way. - dat$oddcol <- factor(dat$col - floor(dat$col / 2) *2 ) - - m3 <- update(m2, yield ~ 1 + oddcol) - m3$loglik # Matches Burgueno table 3, line 3 - - plot(varioGram(m3), xlim=c(0,15), ylim=c(0,9), zlim=c(0,0.5), - main="burgueno.unreplicated - AR1xAR1 + Even/Odd") - # Much better-looking variogram + m3 <- update(m2, yield ~ 1 + oddcol) + m3$loglik # Matches Burgueno table 3, line 3 + + plot(varioGram(m3), xlim=c(0,15), ylim=c(0,9), zlim=c(0,0.5), + main="burgueno.unreplicated - AR1xAR1 + Even/Odd") + # Much better-looking variogram + } } } diff --git a/man/butron.maize.Rd b/man/butron.maize.Rd index 43e1d595..330dfcd1 100644 --- a/man/butron.maize.Rd +++ b/man/butron.maize.Rd @@ -46,7 +46,7 @@ } \source{ - Butron, A and Velasco, P and Ord{\'a}s, A and Malvar, RA (2004). + Butron, A and Velasco, P and Ordas, A and Malvar, RA (2004). Yield evaluation of maize cultivars across environments with different levels of pink stem borer infestation. Crop Science, 44, 741-747. @@ -71,8 +71,8 @@ biplot(princomp(mat), main="butron.maize", cex=.7) # Figure 1 of Butron - if(0){ - # Fixme: This section is broken because synbreed has been removed from CRAN + if(require("asreml", quietly=TRUE)) { + # Here we see if including pedigree information is helpful for a # multi-environment model # Including the pedigree provided little benefit @@ -93,12 +93,7 @@ ped) ped[is.na(ped$male),'male'] <- 0 ped[is.na(ped$female),'female'] <- 0 - - } - - # asreml - if(0){ libs(asreml) ped.ainv <- ainverse(ped) @@ -134,9 +129,10 @@ p0$gen, cex=.5, srt=-45) text(x=min(lims), y=p1par$predicted.value, p1par$gen, cex=.5, col="red") round( cor(p0$predicted.value, p1$predicted.value), 3) + # 0.994 # Including the pedigree provided very little change - } + } } diff --git a/man/cleveland.soil.Rd b/man/cleveland.soil.Rd index a2a1a76a..9269c95c 100644 --- a/man/cleveland.soil.Rd +++ b/man/cleveland.soil.Rd @@ -35,32 +35,32 @@ \examples{ \dontrun{ -library(agridat) -data(cleveland.soil) -dat <- cleveland.soil - -# Similar to Cleveland fig 4.64 -## libs(latticeExtra) -## redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997")) -## levelplot(resistivity ~ easting + northing, data = dat, -## col.regions=redblue, -## panel=panel.levelplot.points, -## aspect=2.4, xlab= "Easting (km)", ylab= "Northing (km)", -## main="cleveland") - -# 2D loess plot. Cleveland fig 4.68 -sg1 <- expand.grid(easting = seq(.15, 1.410, by = .02), - northing = seq(.150, 3.645, by = .02)) -lo1 <- loess(resistivity~easting*northing, data=dat, span = 0.1, degree = 2) -fit1 <- predict(lo1, sg1) -libs(lattice) -redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997")) -levelplot(fit1 ~ sg1$easting * sg1$northing, - col.regions=redblue, - cuts = 9, - aspect=2.4, xlab = "Easting (km)", ylab = "Northing (km)", - main="cleveland.soil - 2D smooth of Resistivity") + library(agridat) + data(cleveland.soil) + dat <- cleveland.soil + # Similar to Cleveland fig 4.64 + ## libs(latticeExtra) + ## redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997")) + ## levelplot(resistivity ~ easting + northing, data = dat, + ## col.regions=redblue, + ## panel=panel.levelplot.points, + ## aspect=2.4, xlab= "Easting (km)", ylab= "Northing (km)", + ## main="cleveland") + + # 2D loess plot. Cleveland fig 4.68 + sg1 <- expand.grid(easting = seq(.15, 1.410, by = .02), + northing = seq(.150, 3.645, by = .02)) + lo1 <- loess(resistivity~easting*northing, data=dat, span = 0.1, degree = 2) + fit1 <- predict(lo1, sg1) + libs(lattice) + redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997")) + levelplot(fit1 ~ sg1$easting * sg1$northing, + col.regions=redblue, + cuts = 9, + aspect=2.4, xlab = "Easting (km)", ylab = "Northing (km)", + main="cleveland.soil - 2D smooth of Resistivity") + # 3D loess plot with data overlaid libs(rgl) bg3d(color = "white") @@ -71,7 +71,7 @@ levelplot(fit1 ~ sg1$easting * sg1$northing, seq(.150, 3.645, by = .02), fit1/100, alpha=0.9, col=rep("wheat", length(fit1)), front="fill", back="fill") - rgl.close() + close3d() } } \keyword{datasets} diff --git a/man/cochran.eelworms.Rd b/man/cochran.eelworms.Rd index 9ac0f77a..119c58cd 100644 --- a/man/cochran.eelworms.Rd +++ b/man/cochran.eelworms.Rd @@ -70,35 +70,35 @@ \examples{ \dontrun{ -library(agridat) -data(cochran.eelworms) -dat <- cochran.eelworms + library(agridat) + data(cochran.eelworms) + dat <- cochran.eelworms -libs(lattice) -splom(dat[ , 5:10], - group=dat$fumigant, auto.key=TRUE, - main="cochran.eelworms") + libs(lattice) + splom(dat[ , 5:10], + group=dat$fumigant, auto.key=TRUE, + main="cochran.eelworms") -libs(desplot) -desplot(dat, fumigant~col*row, text=dose, flip=TRUE, cex=2) + libs(desplot) + desplot(dat, fumigant~col*row, text=dose, flip=TRUE, cex=2) -# Very strong spatial trends -desplot(dat, initial ~ col*row, - flip=TRUE, # aspect unknown - main="cochran.eelworms") + # Very strong spatial trends + desplot(dat, initial ~ col*row, + flip=TRUE, # aspect unknown + main="cochran.eelworms") -# final counts are strongly related to initial counts -libs(lattice) -xyplot(final~initial|factor(dose), data=dat, group=fumigant, - main="cochran.eelworms - by dose (panel) & fumigant", - xlab="Initial worm count", - ylab="Final worm count", auto.key=list(columns=5)) - -# One approach...log transform, use 'initial' as covariate, create 9 treatments -dat <- transform(dat, trt=factor(paste0(fumigant, dose))) -m1 <- aov(log(final) ~ block + trt + log(initial), data=dat) -anova(m1) + # final counts are strongly related to initial counts + libs(lattice) + xyplot(final~initial|factor(dose), data=dat, group=fumigant, + main="cochran.eelworms - by dose (panel) & fumigant", + xlab="Initial worm count", + ylab="Final worm count", auto.key=list(columns=5)) + + # One approach...log transform, use 'initial' as covariate, create 9 treatments + dat <- transform(dat, trt=factor(paste0(fumigant, dose))) + m1 <- aov(log(final) ~ block + trt + log(initial), data=dat) + anova(m1) } } diff --git a/man/cullis.earlygen.Rd b/man/cullis.earlygen.Rd index 8f5ecb1d..f44a8840 100644 --- a/man/cullis.earlygen.Rd +++ b/man/cullis.earlygen.Rd @@ -102,52 +102,53 @@ shapeCross=shape, layers=NULL) dat$mov.avg <- fitted(m0) - libs(asreml) - - # Start with the standard AR1xAR1 analysis - dat <- transform(dat, xf=factor(col), yf=factor(row)) - dat <- dat[order(dat$xf, dat$yf),] - m2 <- asreml(yield ~ weed, data=dat, random= ~gen, - resid = ~ ar1(xf):ar1(yf)) - - # Variogram suggests a polynomial trend - m3 <- update(m2, fixed= yield~weed+pol(col,-1)) - - # Now add a nugget variance - m4 <- update(m3, random= ~ gen + units) - - libs(lucid) - vc(m4) - ## effect component std.error z.ratio bound %ch - ## gen 73780 10420 7.1 P 0 - ## units 30440 8073 3.8 P 0.1 - ## xf:yf(R) 54730 10630 5.1 P 0 - ## xf:yf!xf!cor 0.38 0.115 3.3 U 0 - ## xf:yf!yf!cor 0.84 0.045 19 U 0 - - ## # Predictions from models m3 and m4 are non-estimable. Why? - ## # Use model m2 for predictions - ## predict(m2, classify="gen")$pvals - ## ## gen predicted.value std.error status - ## ## 1 Banks 2723.534 93.14719 Estimable - ## ## 2 Eno008 2981.056 162.85241 Estimable - ## ## 3 Eno009 2978.008 161.57129 Estimable - ## ## 4 Eno010 2821.399 153.96943 Estimable - ## ## 5 Eno011 2991.612 161.53507 Estimable + if(require("asreml", quietly=TRUE)) { + libs(asreml,lucid) + + # Start with the standard AR1xAR1 analysis + dat <- transform(dat, xf=factor(col), yf=factor(row)) + dat <- dat[order(dat$xf, dat$yf),] + m2 <- asreml(yield ~ weed, data=dat, random= ~gen, + resid = ~ ar1(xf):ar1(yf)) + + # Variogram suggests a polynomial trend + m3 <- update(m2, fixed= yield~weed+pol(col,-1)) + + # Now add a nugget variance + m4 <- update(m3, random= ~ gen + units) + + lucid::vc(m4) + ## effect component std.error z.ratio bound %ch + ## gen 73780 10420 7.1 P 0 + ## units 30440 8073 3.8 P 0.1 + ## xf:yf(R) 54730 10630 5.1 P 0 + ## xf:yf!xf!cor 0.38 0.115 3.3 U 0 + ## xf:yf!yf!cor 0.84 0.045 19 U 0 + + ## # Predictions from models m3 and m4 are non-estimable. Why? + ## # Use model m2 for predictions + ## predict(m2, classify="gen")$pvals + ## ## gen predicted.value std.error status + ## ## 1 Banks 2723.534 93.14719 Estimable + ## ## 2 Eno008 2981.056 162.85241 Estimable + ## ## 3 Eno009 2978.008 161.57129 Estimable + ## ## 4 Eno010 2821.399 153.96943 Estimable + ## ## 5 Eno011 2991.612 161.53507 Estimable + + + ## # Compare AR1 with Moving Grid + ## dat$ar1 <- fitted(m2) + ## head(dat[ , c('yield','ar1','mov.avg')]) + ## ## yield ar1 mg + ## ## 1 2652 2467.980 2531.998 + ## ## 11 3394 3071.681 3052.160 + ## ## 21 3148 2826.188 2807.031 + ## ## 31 3426 3026.985 3183.649 + ## ## 41 3555 3070.102 3195.910 + ## ## 51 3453 3006.352 3510.511 + ## pairs(dat[ , c('yield','ar1','mg')]) + } - - ## # Compare AR1 with Moving Grid - ## dat$ar1 <- fitted(m2) - ## head(dat[ , c('yield','ar1','mov.avg')]) - ## ## yield ar1 mg - ## ## 1 2652 2467.980 2531.998 - ## ## 11 3394 3071.681 3052.160 - ## ## 21 3148 2826.188 2807.031 - ## ## 31 3426 3026.985 3183.649 - ## ## 41 3555 3070.102 3195.910 - ## ## 51 3453 3006.352 3510.511 - ## pairs(dat[ , c('yield','ar1','mg')]) - } } \keyword{datasets} diff --git a/man/damesa.maize.Rd b/man/damesa.maize.Rd index 4b7ee0e6..988a1e50 100644 --- a/man/damesa.maize.Rd +++ b/man/damesa.maize.Rd @@ -49,13 +49,13 @@ main="damesa.maize", out1=rep, out2=block, num=gen, cex=1) - if(0){ + if(require("asreml", quietly=TRUE)) { # Fit the single-stage model in Damesa - lib(asreml) + libs(asreml,lucid) m0 <- asreml(data=damesa.maize, - fixed = yield ~ gen, - random = ~ site + gen:site + at(site):rep/block, - residual = ~ dsum( ~ units|site) ) + fixed = yield ~ gen, + random = ~ site + gen:site + at(site):rep/block, + residual = ~ dsum( ~ units|site) ) lucid::vc(m0) # match Damesa table 1 column 3 ## effect component std.error z.ratio bound %ch ## at(site, S1):rep 0.08819 0.1814 0.49 P 0 @@ -73,6 +73,7 @@ ## site_S3!R 1.153 0.2349 4.9 P 0 ## site_S4!R 0.1112 0.03665 3 P 0 } + } } \keyword{datasets} diff --git a/man/digby.jointregression.Rd b/man/digby.jointregression.Rd index 9e0829f8..de453e97 100644 --- a/man/digby.jointregression.Rd +++ b/man/digby.jointregression.Rd @@ -117,9 +117,11 @@ # Using 'mumm' gives similar results, though gen is random and the # coeffecients are shrunk toward 0 a bit. - libs(mumm) - m1 <- mumm(yield ~ -1 + env + mp(gen, env), dat) - round(1 + ranef(m1)$`mp gen:env`,2) + if(require("mumm", quietly=TRUE)) { + libs(mumm) + m1 <- mumm(yield ~ -1 + env + mp(gen, env), dat) + round(1 + ranef(m1)$`mp gen:env`,2) + } } } diff --git a/man/diggle.cow.Rd b/man/diggle.cow.Rd index ec4f8b6d..86c806de 100644 --- a/man/diggle.cow.Rd +++ b/man/diggle.cow.Rd @@ -71,26 +71,28 @@ Bodyweight of cows in a 2-by-2 factorial experiment useOuterStrips(xyplot(weight ~ day|iron*infect, dat, group=animal, type='b', cex=.5, main="diggle.cow")) - + # Scaling dat <- transform(dat, time = (day-122)/10) - libs(asreml) - ## libs(latticeExtra) + if(require("asreml", quietly=TRUE)) { + libs(asreml, latticeExtra) + + ## # Smooth for each animal. No treatment effects. Similar to SAS Output 38.6.9 + + m1 <- asreml(weight ~ 1 + lin(time) + animal + animal:lin(time), data=dat, + random = ~ animal:spl(time)) + p1 <- predict(m1, data=dat, classify="animal:time", + design.points=list(time=seq(0,65.9, length=50))) + p1 <- p1$pvals + p1 <- merge(dat, p1, all=TRUE) # to get iron/infect merged in + foo1 <- xyplot(weight ~ day|iron*infect, dat, group=animal, + main="diggle.cow") + foo2 <- xyplot(predicted.value ~ day|iron*infect, p1, type='l', group=animal) + print(foo1+foo2) + } - ## # Smooth for each animal. No treatment effects. Similar to SAS Output 38.6.9 - - m1 <- asreml(weight ~ 1 + lin(time) + animal + animal:lin(time), data=dat, - random = ~ animal:spl(time)) - p1 <- predict(m1, data=dat, classify="animal:time", - design.points=list(time=seq(0,65.9, length=50))) - p1 <- p1$pvals - p1 <- merge(dat, p1, all=TRUE) # to get iron/infect merged in - foo1 <- xyplot(weight ~ day|iron*infect, dat, group=animal) - foo2 <- xyplot(predicted.value ~ day|iron*infect, p1, type='l', group=animal) - print(foo1+foo2) - } } \keyword{datasets} diff --git a/man/draper.safflower.uniformity.Rd b/man/draper.safflower.uniformity.Rd index cb80a2bd..a89c30bd 100644 --- a/man/draper.safflower.uniformity.Rd +++ b/man/draper.safflower.uniformity.Rd @@ -71,17 +71,16 @@ \examples{ \dontrun{ - library(agridat) - + library(agridat) data(draper.safflower.uniformity) dat4 <- subset(draper.safflower.uniformity, expt=="E4") dat5 <- subset(draper.safflower.uniformity, expt=="E5") - + libs(desplot) desplot(dat4, yield~col*row, flip=TRUE, tick=TRUE, aspect=72/53, # true aspect main="draper.safflower.uniformity (four foot)") - + desplot(dat5, yield~col*row, flip=TRUE, tick=TRUE, aspect=90/46, # true aspect main="draper.safflower.uniformity (five foot)") @@ -91,7 +90,7 @@ dat4 <- subset(dat4, row>1 & row<20) dat4 <- subset(dat4, col>1 & col<17) dat5 <- subset(dat5, row>1 & row<20) - dat5 <- subset(dat5, col<15)\ + dat5 <- subset(dat5, col<15) # Convert gm/plot to pounds/acre. Draper (p. 20) says 1487 pounds/acre mean(dat4$yield) / 453.592 / (3.33*4) * 43560 # 1472 lb/ac diff --git a/man/durban.rowcol.Rd b/man/durban.rowcol.Rd index fcbce7c6..0d3ab2d9 100644 --- a/man/durban.rowcol.Rd +++ b/man/durban.rowcol.Rd @@ -83,32 +83,34 @@ wireframe(p1lo~row+bed, new1, aspect=c(1,.5), main="Field trend") - libs(asreml) + if(require("asreml", quietly=TRUE)) { + libs(asreml) - dat <- transform(dat, rowf=factor(row), bedf=factor(bed)) - dat <- dat[order(dat$rowf, dat$bedf),] + dat <- transform(dat, rowf=factor(row), bedf=factor(bed)) + dat <- dat[order(dat$rowf, dat$bedf),] - m1a1 <- asreml(yield~gen + lin(rowf) + lin(bedf), data=dat, - random=~spl(rowf) + spl(bedf) + units, - family=asr_gaussian(dispersion=1)) - m1a2 <- asreml(yield~gen + lin(rowf) + lin(bedf), data=dat, - random=~spl(rowf) + spl(bedf) + units, - resid = ~ar1(rowf):ar1(bedf)) - m1a2 <- update(m1a2) - m1a3 <- asreml(yield~gen, data=dat, random=~units, - resid = ~ar1(rowf):ar1(bedf)) - - # Figure 7 - libs(lattice) - v7a <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1a3$residuals) - wireframe(gamma ~ x*y, v7a, aspect=c(1,.5)) # Fig 7a - - v7b <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1a2$residuals) - wireframe(gamma ~ x*y, v7b, aspect=c(1,.5)) # Fig 7b + m1a1 <- asreml(yield~gen + lin(rowf) + lin(bedf), data=dat, + random=~spl(rowf) + spl(bedf) + units, + family=asr_gaussian(dispersion=1)) + m1a2 <- asreml(yield~gen + lin(rowf) + lin(bedf), data=dat, + random=~spl(rowf) + spl(bedf) + units, + resid = ~ar1(rowf):ar1(bedf)) + m1a2 <- update(m1a2) + m1a3 <- asreml(yield~gen, data=dat, random=~units, + resid = ~ar1(rowf):ar1(bedf)) + + # Figure 7 + libs(lattice) + v7a <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1a3$residuals) + wireframe(gamma ~ x*y, v7a, aspect=c(1,.5)) # Fig 7a + + v7b <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1a2$residuals) + wireframe(gamma ~ x*y, v7b, aspect=c(1,.5)) # Fig 7b + + v7c <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1lo$residuals) + wireframe(gamma ~ x*y, v7c, aspect=c(1,.5)) # Fig 7c + } - v7c <- asr_varioGram(x=dat$bedf, y=dat$rowf, z=m1lo$residuals) - wireframe(gamma ~ x*y, v7c, aspect=c(1,.5)) # Fig 7c - } } diff --git a/man/durban.splitplot.Rd b/man/durban.splitplot.Rd index c03d44af..77d51353 100644 --- a/man/durban.splitplot.Rd +++ b/man/durban.splitplot.Rd @@ -74,27 +74,28 @@ wireframe(p2lo~row+bed, new2, aspect=c(1,.5), main="durban.splitplot - Field trend") - libs(asreml) + if(require("asreml", quietly=TRUE)) { + libs(asreml,lucid) - # Table 5, variance components. Table 6, F tests - dat <- transform(dat, rowf=factor(row), bedf=factor(bed)) - dat <- dat[order(dat$rowf, dat$bedf),] - m2a2 <- asreml(yield ~ gen*fung, random=~block/fung+units, data=dat, - resid =~ar1v(rowf):ar1(bedf)) - m2a2 <- update(m2a2) - - libs(lucid) - vc(m2a2) - ## effect component std.error z.ratio bound %ch - ## block 0 NA NA B NA - ## block:fung 0.01206 0.01512 0.8 P 0 - ## units 0.02463 0.002465 10 P 0 - ## rowf:bedf(R) 1 NA NA F 0 - ## rowf:bedf!rowf!cor 0.8836 0.03646 24 U 0 - ## rowf:bedf!rowf!var 0.1261 0.04434 2.8 P 0 - ## rowf:bedf!bedf!cor 0.9202 0.02846 32 U 0 - - wald(m2a2) + # Table 5, variance components. Table 6, F tests + dat <- transform(dat, rowf=factor(row), bedf=factor(bed)) + dat <- dat[order(dat$rowf, dat$bedf),] + m2a2 <- asreml(yield ~ gen*fung, random=~block/fung+units, data=dat, + resid =~ar1v(rowf):ar1(bedf)) + m2a2 <- update(m2a2) + + lucid::vc(m2a2) + ## effect component std.error z.ratio bound %ch + ## block 0 NA NA B NA + ## block:fung 0.01206 0.01512 0.8 P 0 + ## units 0.02463 0.002465 10 P 0 + ## rowf:bedf(R) 1 NA NA F 0 + ## rowf:bedf!rowf!cor 0.8836 0.03646 24 U 0 + ## rowf:bedf!rowf!var 0.1261 0.04434 2.8 P 0 + ## rowf:bedf!bedf!cor 0.9202 0.02846 32 U 0 + + wald(m2a2) + } } } diff --git a/man/federer.diagcheck.Rd b/man/federer.diagcheck.Rd index fd46bd1d..b3429aa2 100644 --- a/man/federer.diagcheck.Rd +++ b/man/federer.diagcheck.Rd @@ -41,14 +41,14 @@ \source{ Federer, Walter T. 1998. Recovery of interblock, intergradient, and intervariety information in incomplete block and lattice rectangle - design experiments. \emph{Biometrics}, 54, 471--481. + design experiments. Biometrics, 54, 471--481. https://doi.org/10.2307/3109756 } \references{ Walter T Federer and Russell D Wolfinger, 2003. Augmented Row-Column Design and Trend Analysis, chapter 28 of - \emph{Handbook of Formulas and Software for Plant Geneticists - and Breeders}, Haworth Press. + Handbook of Formulas and Software for Plant Geneticists + and Breeders, Haworth Press. } \examples{ @@ -160,42 +160,44 @@ ## one.12 r2 928.2 30.47 ## one.13 r1 10360 101.8 ## Residual 4127 64.24 - - libs(asreml,lucid) - - m3 <- asreml(yield ~ -1 + trtn, data=dat, - random = ~ r1 + r2 + r4 + r8 + r10 + - c1 + c2 + c3 + c4 + c6 + c8 + - r1:c1 + r1:c2 + r1:c3 + new:gen) - ## coef(m3) - ## # REML cultivar means. Very similar to Federer table 2. - ## rev(sort(round(coef(m3)$fixed[3] + coef(m3)$random[137:256,],0))) - ## ## gen_G060 gen_G021 gen_G011 gen_G099 gen_G002 - ## ## 974 949 945 944 942 - ## ## gen_G118 gen_G058 gen_G035 gen_G111 gen_G120 - ## ## 938 937 937 933 932 - ## ## gen_G046 gen_G061 gen_G082 gen_G038 gen_G090 - ## ## 932 931 927 927 926 - - ## vc(m3) - ## ## effect component std.error z.ratio constr - ## ## r1!r1.var 9201 13720 0.67 pos - ## ## r2!r2.var 241.7 1059 0.23 pos - ## ## r4!r4.var 2269 3915 0.58 pos - ## ## r8!r8.var 1355 2627 0.52 pos - ## ## r10!r10.var 1133 2312 0.49 pos - ## ## c1!c1.var 0.01 0 4.8 bound - ## ## c2!c2.var 5942 8969 0.66 pos - ## ## c3!c3.var 2549 4177 0.61 pos - ## ## c4!c4.var 1792 3106 0.58 pos - ## ## c6!c6.var 1400 2551 0.55 pos - ## ## c8!c8.var 6456 9702 0.67 pos - ## ## r1:c1!r1.var 128000 189700 0.67 pos - ## ## r1:c2!r1.var 58230 90820 0.64 pos - ## ## r1:c3!r1.var 5531 16550 0.33 pos - ## ## new:gen!new.var 2869 1367 2.1 pos - ## ## R!variance 4412 915 4.8 pos + if(require("asreml", quietly=TRUE)) { + libs(asreml,lucid) + + m3 <- asreml(yield ~ -1 + trtn, data=dat, + random = ~ r1 + r2 + r4 + r8 + r10 + + c1 + c2 + c3 + c4 + c6 + c8 + + r1:c1 + r1:c2 + r1:c3 + new:gen) + ## coef(m3) + ## # REML cultivar means. Very similar to Federer table 2. + ## rev(sort(round(coef(m3)$fixed[3] + coef(m3)$random[137:256,],0))) + ## ## gen_G060 gen_G021 gen_G011 gen_G099 gen_G002 + ## ## 974 949 945 944 942 + ## ## gen_G118 gen_G058 gen_G035 gen_G111 gen_G120 + ## ## 938 937 937 933 932 + ## ## gen_G046 gen_G061 gen_G082 gen_G038 gen_G090 + ## ## 932 931 927 927 926 + + ## vc(m3) + ## ## effect component std.error z.ratio constr + ## ## r1!r1.var 9201 13720 0.67 pos + ## ## r2!r2.var 241.7 1059 0.23 pos + ## ## r4!r4.var 2269 3915 0.58 pos + ## ## r8!r8.var 1355 2627 0.52 pos + ## ## r10!r10.var 1133 2312 0.49 pos + ## ## c1!c1.var 0.01 0 4.8 bound + ## ## c2!c2.var 5942 8969 0.66 pos + ## ## c3!c3.var 2549 4177 0.61 pos + ## ## c4!c4.var 1792 3106 0.58 pos + ## ## c6!c6.var 1400 2551 0.55 pos + ## ## c8!c8.var 6456 9702 0.67 pos + ## ## r1:c1!r1.var 128000 189700 0.67 pos + ## ## r1:c2!r1.var 58230 90820 0.64 pos + ## ## r1:c3!r1.var 5531 16550 0.33 pos + ## ## new:gen!new.var 2869 1367 2.1 pos + ## ## R!variance 4412 915 4.8 pos + } + } } \keyword{datasets} diff --git a/man/fisher.barley.Rd b/man/fisher.barley.Rd index 82fc565e..67bbfc8f 100644 --- a/man/fisher.barley.Rd +++ b/man/fisher.barley.Rd @@ -35,17 +35,17 @@ George Fernandez (1991). Analysis of Genotype x Environment Interaction by Stability Estimates. - \emph{Hort Science}, 26, 947-950. + Hort Science, 26, 947-950. F. Yates & W. G. Cochran (1938). The Analysis of Groups of Experiments. - \emph{Journal of Agricultural Science}, 28, 556-580, table 1. + Journal of Agricultural Science, 28, 556-580, table 1. https://doi.org/10.1017/S0021859600050978 G. K. Shukla, 1972. Some statistical aspects of partitioning of genotype-environmental components of variability. - \emph{Heredity}, 29, 237-245. Table 1. + Heredity, 29, 237-245. Table 1. https://doi.org/10.1038/hdy.1972.87 } @@ -97,42 +97,43 @@ ## 4 Trebi 225.52866 ## 5 Velvet 22.73051 + if(require("asreml", quietly=TRUE)) { + # mixed model approach gives similar results (but not identical) + libs(asreml,lucid) - # mixed model approach gives similar results (but not identical) - libs(asreml,lucid) - - dat2 <- dat - dat2 <- dplyr::group_by(dat2, gen,env) - dat2 <- dplyr::summarize(dat2, yield=sum(yield)) # means across years - dat2 <- dplyr::arrange(dat2, gen) + dat2 <- dat + dat2 <- dplyr::group_by(dat2, gen,env) + dat2 <- dplyr::summarize(dat2, yield=sum(yield)) # means across years + dat2 <- dplyr::arrange(dat2, gen) + + # G-side + m1g <- asreml(yield ~ gen, data=dat2, + random = ~ env + at(gen):units, + family=asr_gaussian(dispersion=1.0)) + m1g <- update(m1g) + summary(m1g)$varcomp[-1,1:2]/6 + # component std.error + # at(gen, Manchuria):units 33.8145031 27.22721 + # at(gen, Peatland):units 70.4489092 50.52680 + # at(gen, Svansota):units 25.2728568 21.92919 + # at(gen, Trebi):units 231.6981702 150.80464 + # at(gen, Velvet):units 13.9325646 16.58571 + # units!R 0.1666667 NA + + # R-side estimates = G-side estimate + 0.1666 (resid variance) + m1r <- asreml(yield ~ gen, data=dat2, + random = ~ env, + residual = ~ dsum( ~ units|gen)) + m1r <- update(m1r) + summary(m1r)$varcomp[-1,1:2]/6 + # component std.error + # gen_Manchuria!R 34.00058 27.24871 + # gen_Peatland!R 70.65501 50.58925 + # gen_Svansota!R 25.42022 21.88606 + # gen_Trebi!R 231.85846 150.78756 + # gen_Velvet!R 14.08405 16.55558 + } - # G-side - m1g <- asreml(yield ~ gen, data=dat2, - random = ~ env + at(gen):units, - family=asr_gaussian(dispersion=1.0)) - m1g <- update(m1g) - summary(m1g)$varcomp[-1,1:2]/6 - # component std.error - # at(gen, Manchuria):units 33.8145031 27.22721 - # at(gen, Peatland):units 70.4489092 50.52680 - # at(gen, Svansota):units 25.2728568 21.92919 - # at(gen, Trebi):units 231.6981702 150.80464 - # at(gen, Velvet):units 13.9325646 16.58571 - # units!R 0.1666667 NA - - # R-side estimates = G-side estimate + 0.1666 (resid variance) - m1r <- asreml(yield ~ gen, data=dat2, - random = ~ env, - residual = ~ dsum( ~ units|gen)) - m1r <- update(m1r) - summary(m1r)$varcomp[-1,1:2]/6 - # component std.error - # gen_Manchuria!R 34.00058 27.24871 - # gen_Peatland!R 70.65501 50.58925 - # gen_Svansota!R 25.42022 21.88606 - # gen_Trebi!R 231.85846 150.78756 - # gen_Velvet!R 14.08405 16.55558 - } } \keyword{datasets} diff --git a/man/gartner.corn.Rd b/man/gartner.corn.Rd index b0e00fca..b1ecaed3 100644 --- a/man/gartner.corn.Rd +++ b/man/gartner.corn.Rd @@ -67,85 +67,79 @@ \examples{ \dontrun{ -library(agridat) -data(gartner.corn) -dat <- gartner.corn + library(agridat) + data(gartner.corn) + dat <- gartner.corn -# Calculate yield from mass & moisture -dat <- transform(dat, + # Calculate yield from mass & moisture + dat <- transform(dat, yield=(mass*seconds*(100-moist)/(100-15.5)/56)/(dist*360/6272640)) -# Delete low yield outliers -dat <- subset(dat, yield >50) - -# Group yield into 20 bins for red-gray-blue colors -medy <- median(dat$yield) -ncols <- 20 -wwidth <- 150 -brks <- seq(from = -wwidth/2, to=wwidth/2, length=ncols-1) -brks <- c(-250, brks, 250) # 250 is safe..we cleaned data outside ?(50,450)? -yldbrks <- brks + medy -dat <- transform(dat, yldbin = as.numeric(cut(yield, breaks= yldbrks))) -redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997")) -dat$yieldcolor = redblue(ncols)[dat$yldbin] - -# Polygons for soil map units -# Go to: https://websoilsurvey.nrcs.usda.gov/app/WebSoilSurvey.aspx -# Click: Lat and Long. 43.924, -93.975 -# Click the little AOI rectangle icon. Drag around the field -# In the AOI Properties, enter the Name: Gartner -# Click the tab Soil Map to see map unit symbols, names -# Click: Download Soils Data. Click: Create Download Link. -# Download the zip file and find the soilmu_a_aoi files. - -# Read shape files -libs(sf) -fname <- system.file(package="agridat", "files", "gartner.corn.shp") -shp <- sf::st_read( fname ) - -# Annotate soil map units. Coordinates chosen by hand. -mulabs = data.frame( - name=c("110","319","319","230","105C","110","211","110","211","230","105C"), - x = c(-93.97641, -93.97787, -93.97550, -93.97693, -93.97654, -93.97480, - -93.97375, -93.978284, -93.977617, -93.976715, -93.975929), - y = c(43.92185, 43.92290, 43.92358, 43.92445, 43.92532, 43.92553, - 43.92568, 43.922163, 43.926427, 43.926993, 43.926631) ) -mulabs = st_as_sf( mulabs, coords=c("x","y"), crs=4326) -mulabs = st_transform(mulabs, 2264) - -# Colored points for yield -dat <- st_as_sf(dat, coords=c("long","lat"), crs=4326) - -libs(ggplot2) - -ggplot() + - geom_sf(data=dat, aes(col=yieldcolor) ) + - scale_color_identity() + - geom_sf_label(data=mulabs, aes(label=name), cex=2) + - geom_sf(data=shp["MUSYM"], fill="transparent") + - ggtitle("gartner.corn") + - theme_classic() - - -# Further analyses... -# Trim top and bottom ends of the field -dat <- subset(dat, lat < 43.925850 & lat > 43.921178) -# Identify the soil type for each yield point -# Aggregate points by soil type - - -if(0){ - # Draw a 3D surface. Clearly shows the low drainage area - # Re-run the steps above up, stop before the "Colored points" line. - libs(rgl) - dat <- transform(dat, x=long-min(long), y=lat-min(lat), z=elev-min(elev)) - clear3d() - points3d(dat$x, dat$y, dat$z/50000, - col=redblue(ncols)[dat$yldbin]) - axes3d() - title3d(xlab='x',ylab='y',zlab='elev') - rgl.close() -} + # Delete low yield outliers + dat <- subset(dat, yield >50) + + # Group yield into 20 bins for red-gray-blue colors + medy <- median(dat$yield) + ncols <- 20 + wwidth <- 150 + brks <- seq(from = -wwidth/2, to=wwidth/2, length=ncols-1) + brks <- c(-250, brks, 250) # 250 is safe..we cleaned data outside ?(50,450)? + yldbrks <- brks + medy + dat <- transform(dat, yldbin = as.numeric(cut(yield, breaks= yldbrks))) + redblue <- colorRampPalette(c("firebrick", "lightgray", "#375997")) + dat$yieldcolor = redblue(ncols)[dat$yldbin] + + # Polygons for soil map units + # Go to: https://websoilsurvey.nrcs.usda.gov/app/WebSoilSurvey.aspx + # Click: Lat and Long. 43.924, -93.975 + # Click the little AOI rectangle icon. Drag around the field + # In the AOI Properties, enter the Name: Gartner + # Click the tab Soil Map to see map unit symbols, names + # Click: Download Soils Data. Click: Create Download Link. + # Download the zip file and find the soilmu_a_aoi files. + + # Read shape files + libs(sf) + fname <- system.file(package="agridat", "files", "gartner.corn.shp") + shp <- sf::st_read( fname ) + + # Annotate soil map units. Coordinates chosen by hand. + mulabs = data.frame( + name=c("110","319","319","230","105C","110","211","110","211","230","105C"), + x = c(-93.97641, -93.97787, -93.97550, -93.97693, -93.97654, -93.97480, + -93.97375, -93.978284, -93.977617, -93.976715, -93.975929), + y = c(43.92185, 43.92290, 43.92358, 43.92445, 43.92532, 43.92553, + 43.92568, 43.922163, 43.926427, 43.926993, 43.926631) ) + mulabs = st_as_sf( mulabs, coords=c("x","y"), crs=4326) + mulabs = st_transform(mulabs, 2264) + + # Trim top and bottom ends of the field + dat <- subset(dat, lat < 43.925850 & lat > 43.921178) + # Colored points for yield + dat <- st_as_sf(dat, coords=c("long","lat"), crs=4326) + + libs(ggplot2) + + ggplot() + + geom_sf(data=dat, aes(col=yieldcolor) ) + + scale_color_identity() + + geom_sf_label(data=mulabs, aes(label=name), cex=2) + + geom_sf(data=shp["MUSYM"], fill="transparent") + + ggtitle("gartner.corn") + + theme_classic() + + if(0){ + # Draw a 3D surface. Clearly shows the low drainage area + # Re-run the steps above up, stop before the "Colored points" line. + libs(rgl) + dat <- transform(dat, x=long-min(long), y=lat-min(lat), z=elev-min(elev)) + clear3d() + points3d(dat$x, dat$y, dat$z/50000, + col=redblue(ncols)[dat$yldbin]) + axes3d() + title3d(xlab='x',ylab='y',zlab='elev') + close3d() + } } } diff --git a/man/gilmour.serpentine.Rd b/man/gilmour.serpentine.Rd index 8de832d1..4d8ade31 100644 --- a/man/gilmour.serpentine.Rd +++ b/man/gilmour.serpentine.Rd @@ -65,56 +65,55 @@ # Extreme field trend. Blocking insufficient--needs a spline/smoother # xyplot(yield~col, data=dat, main="gilmour.serpentine") - - # ---------------------------------------------------------------------------- - + if(require("asreml", quietly=TRUE)) { - libs(asreml,lucid) + libs(asreml,lucid) - dat <- transform(dat, rowf=factor(row), colf=factor(10*(col-8))) - dat <- dat[order(dat$rowf, dat$colf), ] # Sort order needed by asreml - - # RCB - m0 <- asreml(yield ~ gen, data=dat, random=~rep) + dat <- transform(dat, rowf=factor(row), colf=factor(10*(col-8))) + dat <- dat[order(dat$rowf, dat$colf), ] # Sort order needed by asreml - # Add AR1 x AR1 - m1 <- asreml(yield ~ gen, data=dat, - resid = ~ar1(rowf):ar1(colf)) - - # Add spline - m2 <- asreml(yield ~ gen + col, data=dat, - random= ~ spl(col) + colf, - resid = ~ar1(rowf):ar1(colf)) - - # Figure 4 shows serpentine spraying - p2 <- predict(m2, data=dat, classify="colf")$pvals - plot(p2$predicted, type='b', xlab="column number", ylab="BLUP") - - # Define column code (due to serpentine spraying) - # Rhelp doesn't like double-percent modulus symbol, so compute by hand - dat <- transform(dat, colcode = factor(dat$col-floor((dat$col-1)/4)*4 -1)) + # RCB + m0 <- asreml(yield ~ gen, data=dat, random=~rep) + + # Add AR1 x AR1 + m1 <- asreml(yield ~ gen, data=dat, + resid = ~ar1(rowf):ar1(colf)) + + # Add spline + m2 <- asreml(yield ~ gen + col, data=dat, + random= ~ spl(col) + colf, + resid = ~ar1(rowf):ar1(colf)) - m3 <- asreml(yield ~ gen + lin(colf) + colcode, data=dat, - random= ~ colf + rowf + spl(colf), - resid = ~ar1(rowf):ar1(colf)) + # Figure 4 shows serpentine spraying + p2 <- predict(m2, data=dat, classify="colf")$pvals + plot(p2$predicted, type='b', xlab="column number", ylab="BLUP") - # Figure 6 shows serpentine row effects - p3 <- predict(m3, data=dat, classify="rowf")$pvals - plot(p3$predicted, type='l', xlab="row number", ylab="BLUP") - text(1:22, p3$predicted, c('L','L','M','R','R','M','L','L', - 'M','R','R','M','L','L','M','R','R','M','L','L','M','R')) + # Define column code (due to serpentine spraying) + # Rhelp doesn't like double-percent modulus symbol, so compute by hand + dat <- transform(dat, colcode = factor(dat$col-floor((dat$col-1)/4)*4 -1)) + + m3 <- asreml(yield ~ gen + lin(colf) + colcode, data=dat, + random= ~ colf + rowf + spl(colf), + resid = ~ar1(rowf):ar1(colf)) - # Define row code (due to serpentine planting). 1=middle, 2=left/right - dat <- transform(dat, rowcode = factor(row)) - levels(dat$rowcode) <- c('2','2','1','2','2','1','2','2','1', - '2','2','1','2','2','1','2','2','1','2','2','1','2') + # Figure 6 shows serpentine row effects + p3 <- predict(m3, data=dat, classify="rowf")$pvals + plot(p3$predicted, type='l', xlab="row number", ylab="BLUP") + text(1:22, p3$predicted, c('L','L','M','R','R','M','L','L', + 'M','R','R','M','L','L','M','R','R','M','L','L','M','R')) + + # Define row code (due to serpentine planting). 1=middle, 2=left/right + dat <- transform(dat, rowcode = factor(row)) + levels(dat$rowcode) <- c('2','2','1','2','2','1','2','2','1', + '2','2','1','2','2','1','2','2','1','2','2','1','2') + + m6 <- asreml(yield ~ gen + lin(colf) + colcode +rowcode, data=dat, + random= ~ colf + rowf + spl(col), + resid = ~ar1(rowf):ar1(colf)) + plot(varioGram(m6), xlim=c(0:17), ylim=c(0,11), zlim=c(0,4000), + main="gilmour.serpentine") + } - m6 <- asreml(yield ~ gen + lin(colf) + colcode +rowcode, data=dat, - random= ~ colf + rowf + spl(col), - resid = ~ar1(rowf):ar1(colf)) - plot(varioGram(m6), xlim=c(0:17), ylim=c(0,11), zlim=c(0,4000), - main="gilmour.serpentine") - } } diff --git a/man/gilmour.slatehall.Rd b/man/gilmour.slatehall.Rd index d0158446..8402580a 100644 --- a/man/gilmour.slatehall.Rd +++ b/man/gilmour.slatehall.Rd @@ -35,8 +35,8 @@ Arthur R Gilmour and Brian R Cullis and Arunas P Verbyla (1997). Accounting for natural and extraneous variation in the analysis of field experiments. - \emph{Journal of Agricultural, Biological, and Environmental - Statistics}, 2, 269-293. + Journal of Agricultural, Biological, and Environmental + Statistics, 2, 269-293. https://doi.org/10.2307/1400446 } @@ -56,29 +56,31 @@ aspect=22.5/40, num=gen, out1=rep, cex=1, main="gilmour.slatehall") - # ---------------------------------------------------------------------------- - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)) { - # Model 4 of Gilmour et al 1997 - dat <- transform(dat, xf=factor(col), yf=factor(row)) - dat <- dat[order(dat$xf, dat$yf), ] - m4 <- asreml(yield ~ gen + lin(row), data=dat, - random = ~ dev(row) + dev(col), - resid = ~ ar1(xf):ar1(yf)) - # coef(m4)$fixed[1] # linear row - # [1] 31.72252 # (sign switch due to row ordering) + libs(asreml,lucid) - vc(m4) - ## effect component std.error z.ratio bound %ch - ## dev(col) 2519 1959 1.3 P 0 - ## dev(row) 20290 10260 2 P 0 - ## xf:yf(R) 23950 4616 5.2 P 0 - ## xf:yf!xf!cor 0.439 0.113 3.9 U 0 - ## xf:yf!yf!cor 0.125 0.117 1.1 U 0 + # Model 4 of Gilmour et al 1997 + dat <- transform(dat, xf=factor(col), yf=factor(row)) + dat <- dat[order(dat$xf, dat$yf), ] + m4 <- asreml(yield ~ gen + lin(row), data=dat, + random = ~ dev(row) + dev(col), + resid = ~ ar1(xf):ar1(yf)) + # coef(m4)$fixed[1] # linear row + # [1] 31.72252 # (sign switch due to row ordering) + + lucid::vc(m4) + ## effect component std.error z.ratio bound %ch + ## dev(col) 2519 1959 1.3 P 0 + ## dev(row) 20290 10260 2 P 0 + ## xf:yf(R) 23950 4616 5.2 P 0 + ## xf:yf!xf!cor 0.439 0.113 3.9 U 0 + ## xf:yf!yf!cor 0.125 0.117 1.1 U 0 + + plot(varioGram(m4), main="gilmour.slatehall") + } - plot(varioGram(m4), main="gilmour.slatehall") - } } \keyword{datasets} diff --git a/man/gomez.multilocsplitplot.Rd b/man/gomez.multilocsplitplot.Rd index 26573f4d..515f9c9f 100644 --- a/man/gomez.multilocsplitplot.Rd +++ b/man/gomez.multilocsplitplot.Rd @@ -31,7 +31,7 @@ } \examples{ - +\dontrun{ library(agridat) data(gomez.multilocsplitplot) dat <- gomez.multilocsplitplot @@ -65,4 +65,4 @@ summary(m1) ## loc:nf:gen 10 1502273 150227 0.3885 } - +} diff --git a/man/hadasch.lettuce.Rd b/man/hadasch.lettuce.Rd index fbb4951e..02ec2d7d 100644 --- a/man/hadasch.lettuce.Rd +++ b/man/hadasch.lettuce.Rd @@ -51,17 +51,17 @@ https://figshare.com/articles/dataset/Lettuce_trial_phenotypic_and_marker_data_/ } \source{ -Hadasch, S., I. Simko, R. J. Hayes, J. O. Ogutu, and H.P. Piepho (2016). -Comparing the predictive abilities of phenotypic and -marker-assisted selection methods in a biparental lettuce population. -Plant Genome 9. -https://doi.org/10.3835/plantgenome2015.03.0014 + Hadasch, S., I. Simko, R. J. Hayes, J. O. Ogutu, and H.P. Piepho (2016). + Comparing the predictive abilities of phenotypic and + marker-assisted selection methods in a biparental lettuce population. + Plant Genome 9. + https://doi.org/10.3835/plantgenome2015.03.0014 } \references{ -Hayes, R. J., Galeano, C. H., Luo, Y., Antonise, R., & Simko, I. (2014). -Inheritance of Decay of Fresh-cut Lettuce in a Recombinant Inbred Line Population from "Salinas 88" × "La Brillante". -J. Amer. Soc. Hort. Sci., 139(4), 388-398. -https://doi.org/10.21273/JASHS.139.4.388 + Hayes, R. J., Galeano, C. H., Luo, Y., Antonise, R., & Simko, I. (2014). + Inheritance of Decay of Fresh-cut Lettuce in a Recombinant Inbred Line Population from "Salinas 88" × "La Brillante". + J. Amer. Soc. Hort. Sci., 139(4), 388-398. + https://doi.org/10.21273/JASHS.139.4.388 } \examples{ \dontrun{ @@ -78,13 +78,15 @@ https://doi.org/10.21273/JASHS.139.4.388 main="hadasch.lettuce") # kinship matrix - tcrossprod(as.matrix(datm[,-1])) %>% head + # head( tcrossprod(as.matrix(datm[,-1])) ) - libs(asreml) - dat <- transform(dat, loc=factor(loc), gen=factor(gen), rep=factor(rep)) - m1 <- asreml(dmr ~ 1 + gen, data=dat, - random = ~ loc + gen:loc + rep:loc) - p1 <- predict(m1, classify="gen")$pvals + if(require("asreml", quietly=TRUE)){ + libs(asreml) + dat <- transform(dat, loc=factor(loc), gen=factor(gen), rep=factor(rep)) + m1 <- asreml(dmr ~ 1 + gen, data=dat, + random = ~ loc + gen:loc + rep:loc) + p1 <- predict(m1, classify="gen")$pvals + } libs(sommer) m2 <- mmer(dmr ~ 0 + gen, data=dat, diff --git a/man/hanks.sprinkler.Rd b/man/hanks.sprinkler.Rd index 051dc380..acee6a28 100644 --- a/man/hanks.sprinkler.Rd +++ b/man/hanks.sprinkler.Rd @@ -41,13 +41,13 @@ Johnson, D. E., Chaudhuri, U. N., and Kanemasu, E. T. (1983). Statistical Analysis of Line-Source Sprinkler Irrigation Experiments and Other Nonrandomized Experiments Using Multivariate Methods. - \emph{Soil Science Society American Journal}, 47, 309-312. + Soil Science Society American Journal, 47, 309-312. Stroup, W. W. (1989). Use of Mixed Model Procedure to Analyze Spatially Correlated Data: An Example Applied to a Line-Source Sprinkler Irrigation Experiment. - \emph{Applications of Mixed Models in Agriculture and Related - Disciplines, Southern Cooperative Series Bulletin No. 343}, 104-122. + Applications of Mixed Models in Agriculture and Related + Disciplines, Southern Cooperative Series Bulletin No. 343, 104-122. SAS Stat User's Guide. https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_mixed_sect038.htm @@ -76,7 +76,6 @@ panel.abline(v=6.5, col='wheat') }) - # ---------------------------------------------------------------------------- ## This is the model from the SAS documentation ## proc mixed; @@ -84,64 +83,34 @@ ## model yield = gen|dir|irr@2; ## random block block*dir block*irr; ## repeated / type=toep(4) sub=block*gen r; - - # ---------------------------------------------------------------------------- - - # asreml 3 - libs(asreml,lucid) - if( utils::packageVersion("asreml") < "4") { - # asreml3 + + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) + dat <- transform(dat, subf=factor(subplot), irrf=factor(irr)) dat <- dat[order(dat$block, dat$gen, dat$subplot),] - m1 <- asreml(yield ~ gen + dir + irrf + gen:dir + gen:irrf + dir:irrf, - data=dat, - random= ~ block + block:dir + block:irrf, - rcov= ~ block:gen:corb(subf, k=4)) - - libs(lucid) - vc(m1) - ## effect component std.error z.ratio constr - ## block!block.var 0.2194 0.2393 0.92 pos - ## block:dir!block.var 0.01768 0.03154 0.56 pos - ## block:irrf!block.var 0.03539 0.03617 0.98 pos - ## R!variance 0.285 0.05086 5.6 pos - ## R!cor1 0.02802 0.1143 0.25 uncon - ## R!cor2 0.005095 0.1278 0.04 uncon - ## R!cor3 -0.3246 0.0905 -3.6 uncon - - ## # convert asreml correlations to SAS covariances - ## round(.2850 * c(1, .02802, .005095, -.3246),4) # res var * (cor1, cor2, cor3) - ## [1] 0.2850 0.0080 0.0015 -0.0925 - } - - # asreml 4 - libs(asreml,lucid) - if( utils::packageVersion("asreml") > "4") { - - dat <- transform(dat, subf=factor(subplot), - irrf=factor(irr)) - dat <- dat[order(dat$block, dat$gen, dat$subplot),] - # In asreml3, we can specify corb(subf, 3) # In asreml4, only corb(subf, 1) runs. corb(subf, 3) says: # Correlation structure is not positive definite m1 <- asreml(yield ~ gen + dir + irrf + gen:dir + gen:irrf + dir:irrf, data=dat, random= ~ block + block:dir + block:irrf, - resid = ~ block:gen:corb(subf, 1)) - + resid = ~ block:gen:corb(subf, 3)) + lucid::vc(m1) - # effect component std.error z.ratio bound %ch - # block 0.194 0.2231 0.87 P 0.5 - # block:dir 0.02729 0.04959 0.55 P 0 - # block:irrf 0.02275 0.0347 0.66 P 0.1 - # block:gen:subf!R 0.3234 0.05921 5.5 P 0 - # block:gen:subf!subf!cor1 0.169 0.09906 1.7 P 0.1 - - } + ## effect component std.error z.ratio bound %ch + ## block 0.2195 0.2378 0.92 P 0.5 + ## block:dir 0.01769 0.03156 0.56 P 0 + ## block:irrf 0.03539 0.0362 0.98 P 0.1 + ## block:gen:subf!R 0.2851 0.05088 5.6 P 0 + ## block:gen:subf!subf!cor1 0.02829 0.1142 0.25 U 0.9 + ## block:gen:subf!subf!cor2 0.004997 0.1278 0.039 U 9.5 + ## block:gen:subf!subf!cor3 -0.3245 0.09044 -3.6 U 0.1 + + } } } diff --git a/man/hansen.multi.uniformity.Rd b/man/hansen.multi.uniformity.Rd index 73e8c4e1..fac9f0b5 100644 --- a/man/hansen.multi.uniformity.Rd +++ b/man/hansen.multi.uniformity.Rd @@ -74,6 +74,7 @@ of plot 10g should probably have been 92 instead of 97. The version of the data in the package uses the changes suggested by Sanders. + Data were typed by K.Wright. } @@ -87,7 +88,7 @@ Eden, T. and E. J. Maskell. (1928). The influence of soil heterogeneity on the growth and yield of successive crops. - Jour of Agricultural Science, 18, 163-185. + Journal of Agricultural Science, 18, 163-185. https://archive.org/stream/in.ernet.dli.2015.25895/2015.25895.Journal-Of-Agricultural-Science-Vol-xviii-1928#page/n175 Sanders, H. G. 1930. @@ -104,25 +105,25 @@ dat <- hansen.multi.uniformity # Field A2: Average across years - libs(dplyr) - dat %>% group_by(row,col) %>% summarize(mn=mean(yield)) %>% dcast(row ~ col, value.var="mn") + libs(dplyr,reshape2) + #dat %>% group_by(row,col) %>% summarize(mn=mean(yield)) %>% dcast(row ~ col, value.var="mn") # Field E2: Match column totals - dat %>% filter(field=="E2") %>% group_by(year,col) %>% summarize(tot=sum(yield)) %>% dcast(year~col, value.var="tot") + #dat %>% filter(field=="E2") %>% group_by(year,col) %>% summarize(tot=sum(yield)) %>% dcast(year~col, value.var="tot") # Heatmaps. Aspect ratio is an educated guess libs(dplyr, desplot) dat <- dat %>% mutate(dat, field=factor(field), year=factor(year)) - dat %>% filter(field=="A2") %>% desplot(yield ~ col*row|year, tick=TRUE, flip=TRUE, aspect=(5*7.4)/(6*7.4)) - dat %>% filter(field=="E2") %>% desplot(yield ~ col*row|year, tick=TRUE, flip=TRUE, aspect=(8*8)/(16*8)) - - # Look at correlation of plots across years - libs(dplyr, lattice) - dat <- dat %>% mutate(plot=paste(row,col)) - filter(dat, field=="A2") %>% acast(plot ~ year, value.var="yield") %>% - splom() - filter(dat, field=="E2") %>% acast(plot ~ year, value.var="yield") %>% - splom() + dat %>% filter(field=="A2") %>% desplot(yield ~ col*row|year, main="hansen.multi.uniformity field A2", tick=TRUE, flip=TRUE, aspect=(5*7.4)/(6*7.4)) + dat %>% filter(field=="E2") %>% desplot(yield ~ col*row|year, main="hansen.multi.uniformity field E2", tick=TRUE, flip=TRUE, aspect=(8*8)/(16*8)) + + # Look at correlation of experimental unit plots across years + libs(dplyr, reshape2, lattice) + dat <- mutate(dat, plot=paste(row,col)) + mat1 <- filter(dat, field=="A2") %>% acast(plot ~ year, value.var="yield") + splom(mat1, main="hansen.multi.uniformity field A2") + mat2 <- filter(dat, field=="E2") %>% acast(plot ~ year, value.var="yield") + splom(mat2, main="hansen.multi.uniformity field A2") } } diff --git a/man/harris.wateruse.Rd b/man/harris.wateruse.Rd index f7148b59..a3a30e88 100644 --- a/man/harris.wateruse.Rd +++ b/man/harris.wateruse.Rd @@ -81,133 +81,129 @@ ## We use pdDiag() to get uncorrelated random effects m1n <- lme(water ~ 1 + ti + ti2, data=d22, na.action=na.omit, random = list(tree=pdDiag(~1+ti+ti2))) - # vc(m1n) + # lucid::vc(m1n) ## effect variance stddev ## (Intercept) 0.2691 0.5188 ## ti 0 0.0000144 ## ti2 0 0.0000039 ## Residual 0.1472 0.3837 - if(0){ - # Various other models with lme4 & asreml - - libs(lme4, lucid) - m1l <- lmer(water ~ 1 + ti + ti2 + (1|tree) + - (0+ti|tree) + (0+ti2|tree), data=d22) - # vc(m1l) - ## grp var1 var2 vcov sdcor - ## tree (Intercept) 0.2691 0.5188 - ## tree.1 ti 0 0 - ## tree.2 ti2 0 0 - ## Residual 0.1472 0.3837 - - - # Once the overall quadratic trend has been removed, there is not - # too much evidence for consecutive observations being correlated - ## d22r <- subset(d22, !is.na(water)) - ## d22r$res <- resid(m1n) - ## xyplot(res ~ day|tree,d22r, - ## as.table=TRUE, type=c('p','smooth'), - ## ylab="residual", - ## main="harris.wateruse - Residuals of individual trees") - ## op <- par(mfrow=c(4,3)) - ## tapply(d22r$res, d22r$tree, acf) - ## par(op) - - # ----- Model 2, add correlation of consecutive measurements - - ## Schabenberger (page 516) adds correlation. - ## Note how the fixed quadratic model is on the "ti = day/100" scale - ## and the correlated observations are on the "day" scale. The - ## only impact this has on the fitted model is to increase the - ## correlation parameter by a factor of 100, which was likely - ## done to get better convergence. - - ## proc mixed data=age2sp2; - ## class tree; - ## model water = ti ti*ti / s ; - ## random intercept /subject=tree s; - ## repeated /subject=tree type=sp(exp)(day); - - ## Same as SAS, use ti for quadratic, day for correlation - m2l <- lme(water ~ 1 + ti + ti2, data=d22, - random = ~ 1|tree, - cor = corExp(form=~ day|tree), - na.action=na.omit) - m2l # Match output 7.27. Same fixef, ranef, variances, exp corr - - # vc(m2l) - ## effect variance stddev - ## (Intercept) 0.2656 0.5154 - ## Residual 0.1541 0.3926 - - # --- - - ## Now use asreml. When I tried rcov=~tree:exp(ti), - ## the estimated parameter value was on the 'boundary', i.e. 0. - ## Changing rcov to the 'day' scale produced a sensible estimate - ## that matched SAS. - ## Note: SAS and asreml use different parameterizations for the correlation - ## SAS uses exp(-d/phi) and asreml uses phi^d. - ## SAS reports 3.79, asreml reports 0.77, and exp(-1/3.7945) = 0.7683274 - ## Note: normally a quadratic would be included as 'pol(day,2)' + # Various other models with lme4 & asreml + + libs(lme4, lucid) + m1l <- lmer(water ~ 1 + ti + ti2 + (1|tree) + + (0+ti|tree) + (0+ti2|tree), data=d22) + # lucid::vc(m1l) + ## grp var1 var2 vcov sdcor + ## tree (Intercept) 0.2691 0.5188 + ## tree.1 ti 0 0 + ## tree.2 ti2 0 0 + ## Residual 0.1472 0.3837 + + + # Once the overall quadratic trend has been removed, there is not + # too much evidence for consecutive observations being correlated + ## d22r <- subset(d22, !is.na(water)) + ## d22r$res <- resid(m1n) + ## xyplot(res ~ day|tree,d22r, + ## as.table=TRUE, type=c('p','smooth'), + ## ylab="residual", + ## main="harris.wateruse - Residuals of individual trees") + ## op <- par(mfrow=c(4,3)) + ## tapply(d22r$res, d22r$tree, acf) + ## par(op) + + # ----- Model 2, add correlation of consecutive measurements + + ## Schabenberger (page 516) adds correlation. + ## Note how the fixed quadratic model is on the "ti = day/100" scale + ## and the correlated observations are on the "day" scale. The + ## only impact this has on the fitted model is to increase the + ## correlation parameter by a factor of 100, which was likely + ## done to get better convergence. + + ## proc mixed data=age2sp2; + ## class tree; + ## model water = ti ti*ti / s ; + ## random intercept /subject=tree s; + ## repeated /subject=tree type=sp(exp)(day); + + ## Same as SAS, use ti for quadratic, day for correlation + m2l <- lme(water ~ 1 + ti + ti2, data=d22, + random = ~ 1|tree, + cor = corExp(form=~ day|tree), + na.action=na.omit) + m2l # Match output 7.27. Same fixef, ranef, variances, exp corr + + # lucid::vc(m2l) + ## effect variance stddev + ## (Intercept) 0.2656 0.5154 + ## Residual 0.1541 0.3926 + + # --- + + ## Now use asreml. When I tried rcov=~tree:exp(ti), + ## the estimated parameter value was on the 'boundary', i.e. 0. + ## Changing rcov to the 'day' scale produced a sensible estimate + ## that matched SAS. + ## Note: SAS and asreml use different parameterizations for the correlation + ## SAS uses exp(-d/phi) and asreml uses phi^d. + ## SAS reports 3.79, asreml reports 0.77, and exp(-1/3.7945) = 0.7683274 + ## Note: normally a quadratic would be included as 'pol(day,2)' + + if(require("asreml", quietly=TRUE)){ libs(asreml) d22 <- d22[order(d22$tree, d22$day),] m2a <- asreml(water ~ 1 + ti + ti2, data=d22, random = ~ tree, residual=~tree:exp(day)) - vc(m2a) + lucid::vc(m2a) ## effect component std.error z.ratio constr ## tree!tree.var 0.2656 0.1301 2 pos ## R!variance 0.1541 0.01611 9.6 pos ## R!day.pow 0.7683 0.04191 18 uncon + } + + # ----- Model 3. Full model for all species/ages. Schabenberger p. 518 + + ## /* Continuous AR(1) autocorrelations included */ + ## proc mixed data=wateruse; + ## class age species tree; + ## model water = age*species age*species*ti age*species*ti*ti / noint s; + ## random intercept ti / subject=age*species*tree s; + ## repeated / subject=age*species*tree type=sp(exp)(day); + + + m3l <- lme(water ~ 0 + age:species + age:species:ti + age:species:ti2, + data=dat, na.action=na.omit, + random = list(tree=pdDiag(~1+ti)), + cor = corExp(form=~ day|tree) ) + + m3l # Match Schabenberger output 7.27. Same fixef, ranef, variances, exp corr + # lucid::vc(m3l) + ## effect variance stddev + ## (Intercept) 0.1549 0.3936 + ## ti 0.02785 0.1669 + ## Residual 0.16 0.4 - # ----- Model 3. Full model for all species/ages. Schabenberger p. 518 - - ## /* Continuous AR(1) autocorrelations included */ - ## proc mixed data=wateruse; - ## class age species tree; - ## model water = age*species age*species*ti age*species*ti*ti / noint s; - ## random intercept ti / subject=age*species*tree s; - ## repeated / subject=age*species*tree type=sp(exp)(day); - - - m3l <- lme(water ~ 0 + age:species + age:species:ti + age:species:ti2, - data=dat, na.action=na.omit, - random = list(tree=pdDiag(~1+ti)), - cor = corExp(form=~ day|tree) ) - - m3l # Match Schabenberger output 7.27. Same fixef, ranef, variances, exp corr - - # vc(m3l) - ## effect variance stddev - ## (Intercept) 0.1549 0.3936 - ## ti 0.02785 0.1669 - ## Residual 0.16 0.4 - - # --- asreml + # --- asreml + if(require("asreml", quietly=TRUE)){ dat <- dat[order(dat$tree,dat$day),] m3a <- asreml(water ~ 0 + age:species + age:species:ti + age:species:ti2, data=dat, random = ~ age:species:tree + age:species:tree:ti, residual = ~ tree:exp(day) ) - # code for asremlr 4 - ## m3a <- asreml(water ~ 0 + age:species + age:species:ti + age:species:ti2, - ## data=dat, - ## random = ~ age:species:tree + age:species:tree:ti, - ## resid = ~ tree:exp(day) ) - # vc(m3a) # Note: day.pow = .8091 = exp(-1/4.7217) + # lucid::vc(m3a) # Note: day.pow = .8091 = exp(-1/4.7217) ## effect component std.error z.ratio constr ## age:species:tree!age.var 0.1549 0.07192 2.2 pos ## age:species:tree:ti!age.var 0.02785 0.01343 2.1 pos ## R!variance 0.16 0.008917 18 pos ## R!day.pow 0.8091 0.01581 51 uncon - } } diff --git a/man/hayman.tobacco.Rd b/man/hayman.tobacco.Rd index 1f3080cc..fce28b3c 100644 --- a/man/hayman.tobacco.Rd +++ b/man/hayman.tobacco.Rd @@ -41,7 +41,7 @@ \source{ B. I. Hayman (1954a). The Analysis of Variance of Diallel Tables. - \emph{Biometrics}, 10, 235-244. Table 5, page 241. + Biometrics, 10, 235-244. Table 5, page 241. https://doi.org/10.2307/3001877 Hayman, B.I. (1954b). The theory and analysis of diallel crosses. @@ -56,13 +56,13 @@ Mohring, Melchinger, Piepho. (2011). REML-Based Diallel Analysis. - \emph{Crop Science}, 51, 470-478. + Crop Science, 51, 470-478. # For 1952 data C. Clark Cockerham and B. S. Weir. (1977). Quadratic analyses of reciprocal crosses. - \emph{Biometrics} 33, 187-203. Appendix C. + Biometrics, 33, 187-203. Appendix C. Andrea Onofri, Niccolo Terzaroli, Luigi Russi (2020). Linear models for diallel crosses: A review with R functions. diff --git a/man/heady.fertilizer.Rd b/man/heady.fertilizer.Rd index ddb56a12..60119dc5 100644 --- a/man/heady.fertilizer.Rd +++ b/man/heady.fertilizer.Rd @@ -99,7 +99,7 @@ summary(m4) with(d2, plot3d(K,P,yield)) with(d3, plot3d(K,P,yield)) with(d4, plot3d(N,P,yield)) # Mostly linear in both N and P - rgl.close() + close3d() } } diff --git a/man/hildebrand.systems.Rd b/man/hildebrand.systems.Rd index 614a8419..53f13372 100644 --- a/man/hildebrand.systems.Rd +++ b/man/hildebrand.systems.Rd @@ -72,46 +72,48 @@ # ---------- - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) - # Environmental variance model, unstructured correlations + # Environmental variance model, unstructured correlations - dat <- dat[order(dat$system, dat$farm),] - m1 <- asreml(yield ~ system, data=dat, - resid = ~us(system):farm) - - # Means, table 5 - ## predict(m1, data=dat, classify="system")$pvals - ## system pred.value std.error est.stat - ## CCA 1.164 0.2816 Estimable - ## CCAF 2.657 0.3747 Estimable - ## LM 1.35 0.1463 Estimable - ## LMF 2.7 0.2561 Estimable - - # Variances, table 5 - # vc(m1)[c(2,4,7,11),] - ## effect component std.error z.ratio constr - ## R!system.CCA:CCA 1.11 0.4354 2.5 pos - ## R!system.CCAF:CCAF 1.966 0.771 2.5 pos - ## R!system.LM:LM 0.2996 0.1175 2.5 pos - ## R!system.LMF:LMF 0.9185 0.3603 2.5 pos - - # Stability variance model - m2 <- asreml(yield ~ system, data=dat, - random = ~ farm, - resid = ~ dsum( ~ units|system)) - m2 <- update(m2) - # predict(m2, data=dat, classify="system")$pvals + dat <- dat[order(dat$system, dat$farm),] + m1 <- asreml(yield ~ system, data=dat, + resid = ~us(system):farm) + + # Means, table 5 + ## predict(m1, data=dat, classify="system")$pvals + ## system pred.value std.error est.stat + ## CCA 1.164 0.2816 Estimable + ## CCAF 2.657 0.3747 Estimable + ## LM 1.35 0.1463 Estimable + ## LMF 2.7 0.2561 Estimable + + # Variances, table 5 + # lucid::vc(m1)[c(2,4,7,11),] + ## effect component std.error z.ratio constr + ## R!system.CCA:CCA 1.11 0.4354 2.5 pos + ## R!system.CCAF:CCAF 1.966 0.771 2.5 pos + ## R!system.LM:LM 0.2996 0.1175 2.5 pos + ## R!system.LMF:LMF 0.9185 0.3603 2.5 pos + + # Stability variance model + m2 <- asreml(yield ~ system, data=dat, + random = ~ farm, + resid = ~ dsum( ~ units|system)) + m2 <- update(m2) + # predict(m2, data=dat, classify="system")$pvals + + ## # Variances, table 6 + # lucid::vc(m2) + ## effect component std.error z.ratio bound %ch + ## farm 0.2998 0.1187 2.5 P 0 + ## system_CCA!R 0.4133 0.1699 2.4 P 0 + ## system_CCAF!R 1.265 0.5152 2.5 P 0 + ## system_LM!R 0.0003805 0.05538 0.0069 P 1.5 + ## system_LMF!R 0.5294 0.2295 2.3 P 0 + } - ## # Variances, table 6 - # vc(m2) - ## effect component std.error z.ratio bound %ch - ## farm 0.2998 0.1187 2.5 P 0 - ## system_CCA!R 0.4133 0.1699 2.4 P 0 - ## system_CCAF!R 1.265 0.5152 2.5 P 0 - ## system_LM!R 0.0003805 0.05538 0.0069 P 1.5 - ## system_LMF!R 0.5294 0.2295 2.3 P 0 - } } diff --git a/man/ilri.sheep.Rd b/man/ilri.sheep.Rd index 07ee84d0..d2816268 100644 --- a/man/ilri.sheep.Rd +++ b/man/ilri.sheep.Rd @@ -121,27 +121,29 @@ # ---------- - libs(asreml,lucid) - # case study page 4.20 - m1 <- asreml(weanwt ~ year + sex + weanage + dl + dq + ramgen + ewegen, - data=dat) - # wald(m1) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) + # case study page 4.20 + m1 <- asreml(weanwt ~ year + sex + weanage + dl + dq + ramgen + ewegen, + data=dat) + # wald(m1) - # case study page 4.26 - m2 <- asreml(weanwt ~ year + sex + weanage + dl + dq + ramgen + ewegen, - random = ~ ram + ewe, data=dat) - # wald(m2) - - # case study page 4.37, year means - # predict(m2, data=dat, classify="year") - ## year predicted.value standard.error est.status - ## 1 91 12.638564 0.2363652 Estimable - ## 2 92 11.067659 0.2285252 Estimable - ## 3 93 11.561932 0.1809891 Estimable - ## 4 94 9.636058 0.2505478 Estimable - ## 5 95 9.350247 0.2346849 Estimable - ## 6 96 10.188482 0.2755387 Estimable + # case study page 4.26 + m2 <- asreml(weanwt ~ year + sex + weanage + dl + dq + ramgen + ewegen, + random = ~ ram + ewe, data=dat) + # wald(m2) + # case study page 4.37, year means + # predict(m2, data=dat, classify="year") + ## year predicted.value standard.error est.status + ## 1 91 12.638564 0.2363652 Estimable + ## 2 92 11.067659 0.2285252 Estimable + ## 3 93 11.561932 0.1809891 Estimable + ## 4 94 9.636058 0.2505478 Estimable + ## 5 95 9.350247 0.2346849 Estimable + ## 6 96 10.188482 0.2755387 Estimable + } + } } \keyword{datasets} diff --git a/man/john.alpha.Rd b/man/john.alpha.Rd index d2a5d62e..eb811b74 100644 --- a/man/john.alpha.Rd +++ b/man/john.alpha.Rd @@ -101,40 +101,39 @@ # ---------- - # asreml - - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) # Heritability calculation of Piepho & Mohring, Example 1 - - m3 <- asreml(yield ~ 1 + rep, data=dat, random=~ gen + rep:block) - sg2 <- summary(m3)$varcomp['gen','component'] # .142902 - - # Average variance of a difference of two adjusted means (BLUP) - - p3 <- predict(m3, data=dat, classify="gen", sed=TRUE) - # Matrix of pair-wise SED values, squared - vdiff <- p3$sed^2 - # Average variance of two DIFFERENT means (using lower triangular of vdiff) - vblup <- mean(vdiff[lower.tri(vdiff)]) # .05455038 - - # Note that without sed=TRUE, asreml reports square root of the average variance - # of a difference between the variety means, so the following gives the same value - # predict(m3, data=dat, classify="gen")$avsed ^ 2 # .05455038 - - # Average variance of a difference of two adjusted means (BLUE) - m4 <- asreml(yield ~ 1 + gen + rep, data=dat, random = ~ rep:block) - p4 <- predict(m4, data=dat, classify="gen", sed=TRUE) - vdiff <- p4$sed^2 - vblue <- mean(vdiff[lower.tri(vdiff)]) # .07010875 - # Again, could use predict(m4, data=dat, classify="gen")$avsed ^ 2 - - # H^2 Ad-hoc measure of heritability - sg2 / (sg2 + vblue/2) # .803 - - # H^2c Similar measure proposed by Cullis. - 1-(vblup / (2*sg2)) # .809 - + + m3 <- asreml(yield ~ 1 + rep, data=dat, random=~ gen + rep:block) + sg2 <- summary(m3)$varcomp['gen','component'] # .142902 + + # Average variance of a difference of two adjusted means (BLUP) + + p3 <- predict(m3, data=dat, classify="gen", sed=TRUE) + # Matrix of pair-wise SED values, squared + vdiff <- p3$sed^2 + # Average variance of two DIFFERENT means (using lower triangular of vdiff) + vblup <- mean(vdiff[lower.tri(vdiff)]) # .05455038 + + # Note that without sed=TRUE, asreml reports square root of the average variance + # of a difference between the variety means, so the following gives the same value + # predict(m3, data=dat, classify="gen")$avsed ^ 2 # .05455038 + + # Average variance of a difference of two adjusted means (BLUE) + m4 <- asreml(yield ~ 1 + gen + rep, data=dat, random = ~ rep:block) + p4 <- predict(m4, data=dat, classify="gen", sed=TRUE) + vdiff <- p4$sed^2 + vblue <- mean(vdiff[lower.tri(vdiff)]) # .07010875 + # Again, could use predict(m4, data=dat, classify="gen")$avsed ^ 2 + + # H^2 Ad-hoc measure of heritability + sg2 / (sg2 + vblue/2) # .803 + + # H^2c Similar measure proposed by Cullis. + 1-(vblup / (2*sg2)) # .809 + } # ---------- # lme4 to calculate Cullis H2 diff --git a/man/kalamkar.wheat.uniformity.Rd b/man/kalamkar.wheat.uniformity.Rd index 6a575ae4..8bcf2825 100644 --- a/man/kalamkar.wheat.uniformity.Rd +++ b/man/kalamkar.wheat.uniformity.Rd @@ -79,19 +79,21 @@ # ---------- - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) - # Show the negative correlation between rows - dat <- transform(dat, - rowf=factor(row), colf=factor(col)) - dat <- dat[order(dat$rowf, dat$colf),] - m1 = asreml(yield ~ 1, data=dat, resid= ~ ar1(rowf):ar1(colf)) - vc(m1) - ## effect component std.error z.ratio bound pctch - ## rowf:colf!R 81.53 3.525 23 P 0 - ## rowf:colf!rowf!cor -0.09464 0.0277 -3.4 U 0.1 - ## rowf:colf!colf!cor 0.2976 0.02629 11 U 0.1 - + # Show the negative correlation between rows + dat <- transform(dat, + rowf=factor(row), colf=factor(col)) + dat <- dat[order(dat$rowf, dat$colf),] + m1 = asreml(yield ~ 1, data=dat, resid= ~ ar1(rowf):ar1(colf)) + lucid::vc(m1) + ## effect component std.error z.ratio bound pctch + ## rowf:colf!R 81.53 3.525 23 P 0 + ## rowf:colf!rowf!cor -0.09464 0.0277 -3.4 U 0.1 + ## rowf:colf!colf!cor 0.2976 0.02629 11 U 0.1 + } + } } \keyword{datasets} diff --git a/man/kempton.barley.uniformity.Rd b/man/kempton.barley.uniformity.Rd index ea30c50b..127e4e02 100644 --- a/man/kempton.barley.uniformity.Rd +++ b/man/kempton.barley.uniformity.Rd @@ -60,20 +60,22 @@ # ---------- - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) - dat <- transform(dat, xf = factor(col), yf=factor(row)) - m1 <- asreml(yield ~ 1, data=dat, resid = ~ar1(xf):ar1(yf)) + dat <- transform(dat, xf = factor(col), yf=factor(row)) + m1 <- asreml(yield ~ 1, data=dat, resid = ~ar1(xf):ar1(yf)) - # vc(m1) - ## effect component std.error z.ratio bound %ch - ## xf:yf!R 0.1044 0.02197 4.7 P 0 - ## xf:yf!xf!cor 0.2458 0.07484 3.3 U 0 - ## xf:yf!yf!cor 0.8186 0.03821 21 U 0 + # lucid::vc(m1) + ## effect component std.error z.ratio bound %ch + ## xf:yf!R 0.1044 0.02197 4.7 P 0 + ## xf:yf!xf!cor 0.2458 0.07484 3.3 U 0 + ## xf:yf!yf!cor 0.8186 0.03821 21 U 0 + + # asreml estimates auto-regression correlations of 0.25, 0.82 + # Kempton estimated auto-regression coefficients b1=0.10, b2=0.91 + } - # asreml estimates auto-regression correlations of 0.25, 0.82 - # Kempton estimated auto-regression coefficients b1=0.10, b2=0.91 - # ---------- if(0){ diff --git a/man/kempton.slatehall.Rd b/man/kempton.slatehall.Rd index 5d48d081..2fed57c0 100644 --- a/man/kempton.slatehall.Rd +++ b/man/kempton.slatehall.Rd @@ -75,25 +75,24 @@ # ---------- - # asreml3 & asreml4 - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) - # Incomplete block model of Gilmour et al 1995 - dat <- transform(dat, xf=factor(col), yf=factor(row)) - m2 <- asreml(yield ~ gen, random = ~ rep/(xf+yf), data=dat) - - vc(m2) - ## effect component std.error z.ratio constr - ## rep!rep.var 4262 6890 0.62 pos - ## rep:xf!rep.var 14810 4865 3 pos - ## rep:yf!rep.var 15600 5091 3.1 pos - ## R!variance 8062 1340 6 pos + # Incomplete block model of Gilmour et al 1995 + dat <- transform(dat, xf=factor(col), yf=factor(row)) + m2 <- asreml(yield ~ gen, random = ~ rep/(xf+yf), data=dat) + + lucid::vc(m2) + ## effect component std.error z.ratio constr + ## rep!rep.var 4262 6890 0.62 pos + ## rep:xf!rep.var 14810 4865 3 pos + ## rep:yf!rep.var 15600 5091 3.1 pos + ## R!variance 8062 1340 6 pos + + # Table 4 + # asreml4 + # predict(m2, data=dat, classify="gen")$pvals + } - # Table 4 - # asreml3 - # predict(m2, data=dat, classify="gen")$predictions$pvals - # asreml4 - # predict(m2, data=dat, classify="gen")$pvals - } } diff --git a/man/kenward.cattle.Rd b/man/kenward.cattle.Rd index 0373bf83..bf4a567b 100644 --- a/man/kenward.cattle.Rd +++ b/man/kenward.cattle.Rd @@ -99,39 +99,41 @@ # ASREML approach to model. Not final by any means. # Maybe a spline curve for each treatment, plus random deviations for each time - libs(asreml) - m1 <- asreml(weight ~ 1 + lin(day) + # overall line - trt + trt:lin(day), # different line for each treatment - data=dat, - random = ~ spl(day) + # overall spline - trt:spl(day) + # different spline for each treatment - dev(day) + trt:dev(day) ) # non-spline deviation at each time*trt + if(require("asreml", quietly=TRUE)){ + libs(asreml) + m1 <- asreml(weight ~ 1 + lin(day) + # overall line + trt + trt:lin(day), # different line for each treatment + data=dat, + random = ~ spl(day) + # overall spline + trt:spl(day) + # different spline for each treatment + dev(day) + trt:dev(day) ) # non-spline deviation at each time*trt + + p1 <- predict(m1, data=dat, classify="trt:day") + p1 <- p1$pvals + + foo2 <- xyplot(predicted.value ~ day|trt, p1, type='l', lwd=2, lty=1, col="black") - p1 <- predict(m1, data=dat, classify="trt:day") - p1 <- p1$pvals - - foo2 <- xyplot(predicted.value ~ day|trt, p1, type='l', lwd=2, lty=1, col="black") - - libs(latticeExtra) - print(foo1 + foo2) + libs(latticeExtra) + print(foo1 + foo2) - # Not much evidence for treatment differences + # Not much evidence for treatment differences - # wald(m1) - ## Df Sum of Sq Wald statistic Pr(Chisq) - ## (Intercept) 1 37128459 139060 <2e-16 *** - ## trt 1 455 2 0.1917 - ## lin(day) 1 570798 2138 <2e-16 *** - ## trt:lin(day) 1 283 1 0.3031 - ## residual (MS) 267 + # wald(m1) + ## Df Sum of Sq Wald statistic Pr(Chisq) + ## (Intercept) 1 37128459 139060 <2e-16 *** + ## trt 1 455 2 0.1917 + ## lin(day) 1 570798 2138 <2e-16 *** + ## trt:lin(day) 1 283 1 0.3031 + ## residual (MS) 267 - # vc(m1) - ## effect component std.error z.ratio constr - ## spl(day) 25.29 24.09 1 pos - ## dev(day) 1.902 4.923 0.39 pos - ## trt:spl(day)!trt.var 0.00003 0.000002 18 bnd - ## trt:dev(day)!trt.var 0.00003 0.000002 18 bnd - ## R!variance 267 14.84 18 pos + # lucid::vc(m1) + ## effect component std.error z.ratio constr + ## spl(day) 25.29 24.09 1 pos + ## dev(day) 1.902 4.923 0.39 pos + ## trt:spl(day)!trt.var 0.00003 0.000002 18 bnd + ## trt:dev(day)!trt.var 0.00003 0.000002 18 bnd + ## R!variance 267 14.84 18 pos + } } } diff --git a/man/lonnquist.maize.Rd b/man/lonnquist.maize.Rd index db4821e7..3fd86fce 100644 --- a/man/lonnquist.maize.Rd +++ b/man/lonnquist.maize.Rd @@ -93,24 +93,25 @@ # ---------- - # asreml - # Mohring 2011 used 6 varieties to calculate GCA & SCA - # Matches Table 3, column 2 - d2 <- subset(dat, is.element(p1, c("M","H","G","B","K","K2")) & - is.element(p2, c("M","H","G","B","K","K2"))) - d2 <- droplevels(d2) - libs(asreml,lucid) - m2 <- asreml(yield~ 1, data=d2, random = ~ p1 + and(p2)) - vc(m2) - ## effect component std.error z.ratio con - ## p1!p1.var 3.865 3.774 1 Positive - ## R!variance 15.93 5.817 2.7 Positive + if(require("asreml", quietly=TRUE)){ + # Mohring 2011 used 6 varieties to calculate GCA & SCA + # Matches Table 3, column 2 + d2 <- subset(dat, is.element(p1, c("M","H","G","B","K","K2")) & + is.element(p2, c("M","H","G","B","K","K2"))) + d2 <- droplevels(d2) + libs(asreml,lucid) + m2 <- asreml(yield~ 1, data=d2, random = ~ p1 + and(p2)) + lucid::vc(m2) + ## effect component std.error z.ratio con + ## p1!p1.var 3.865 3.774 1 Positive + ## R!variance 15.93 5.817 2.7 Positive + + # Calculate GCA effects + m3 <- asreml(yield~ p1 + and(p2), data=d2) + coef(m3)$fixed-1.462 + # Matches Gardner 1966, Table 5, Griffing method + } - # Calculate GCA effects - m3 <- asreml(yield~ p1 + and(p2), data=d2) - coef(m3)$fixed-1.462 - # Matches Gardner 1966, Table 5, Griffing method - } } \keyword{datasets} diff --git a/man/payne.wheat.Rd b/man/payne.wheat.Rd index c7891c3f..0d173aa9 100644 --- a/man/payne.wheat.Rd +++ b/man/payne.wheat.Rd @@ -50,53 +50,54 @@ \examples{ \dontrun{ -library(agridat) -data(payne.wheat) -dat <- payne.wheat + library(agridat) + data(payne.wheat) + dat <- payne.wheat -libs(asreml) - -# make factors -dat <- transform(dat, - rotf = factor(rotation), - yrf = factor(year), - nitrof = factor(nitro)) - -# visualize the response to nitrogen -libs(lattice) -# Why does Payne use nitrogen factor, when it is an obvious polynomial term? -# Probably doesn't matter too much. -xyplot(yield ~ nitro|yrf, dat, - groups=rotf, type='b', - auto.key=list(columns=6), - main="payne.wheat") - -# What are the long-term trends? Yields are decreasing -xyplot(yield ~ year | rotf, data=dat, groups=nitrof, - type='l', auto.key=list(columns=4)) + # make factors + dat <- transform(dat, + rotf = factor(rotation), + yrf = factor(year), + nitrof = factor(nitro)) + + # visualize the response to nitrogen + libs(lattice) + # Why does Payne use nitrogen factor, when it is an obvious polynomial term? + # Probably doesn't matter too much. + xyplot(yield ~ nitro|yrf, dat, + groups=rotf, type='b', + auto.key=list(columns=6), + main="payne.wheat") + + # What are the long-term trends? Yields are decreasing + xyplot(yield ~ year | rotf, data=dat, groups=nitrof, + type='l', auto.key=list(columns=4)) + if(require("asreml", quietly=TRUE)){ + libs(asreml) # Model 5: drop 3-way interaction and return to pol function (easier prediction) -m5 <- asreml(yield ~ rotf * nitrof * pol(year,2) - - (rotf:nitrof:pol(year,2)), - data=dat, - random = ~yrf, - residual = ~ dsum( ~ units|yrf)) -summary(m5)$varcomp # Table 7 of Payne -# lucid::vc(m5) - -# Table 8 of Payne -wald(m5, denDF="default") - -# Predictions of three-way interactions from final model -p5 <- predict(m5, classify="rotf:nitrof:year") -p5 <- p5$pvals # Matches Payne table 8 -head(p5) - -# Plot the predictions. Matches Payne figure 1 -xyplot(predicted.value ~ year | rotf, data=p5, - groups=nitrof, - ylab="yield t/ha", type='l', auto.key=list(columns=5)) + m5 <- asreml(yield ~ rotf * nitrof * pol(year,2) - + (rotf:nitrof:pol(year,2)), + data=dat, + random = ~yrf, + residual = ~ dsum( ~ units|yrf)) + summary(m5)$varcomp # Table 7 of Payne + # lucid::vc(m5) + # Table 8 of Payne + wald(m5, denDF="default") + + # Predictions of three-way interactions from final model + p5 <- predict(m5, classify="rotf:nitrof:year") + p5 <- p5$pvals # Matches Payne table 8 + head(p5) + + # Plot the predictions. Matches Payne figure 1 + xyplot(predicted.value ~ year | rotf, data=p5, + groups=nitrof, + ylab="yield t/ha", type='l', auto.key=list(columns=5)) + } + } } \keyword{datasets} diff --git a/man/perry.springwheat.Rd b/man/perry.springwheat.Rd index 88e83f34..7bbb609c 100644 --- a/man/perry.springwheat.Rd +++ b/man/perry.springwheat.Rd @@ -98,20 +98,22 @@ ## (Intercept) yorc ## 5.87492444 0.05494464 - libs(asreml,lucid) - m3a <- asreml(y ~ 1 + yorc, data=dat, random = ~ env + env:yorc) - vc(m3) - ## grp var1 var2 vcov sdcor - ## env (Intercept) 11.61 3.407 - ## env.1 yorc 0.00063 0.02511 - ## Residual 3.551 1.884 - - vc(m3a) - ## effect component std.error z.ratio con - ## env!env.var 11.61 4.385 2.6 Positive - ## env:yorc!env.var 0.00063 0.000236 2.7 Positive - ## R!variance 3.551 0.2231 16 Positive - + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) + m3a <- asreml(y ~ 1 + yorc, data=dat, random = ~ env + env:yorc) + lucid::vc(m3) + ## grp var1 var2 vcov sdcor + ## env (Intercept) 11.61 3.407 + ## env.1 yorc 0.00063 0.02511 + ## Residual 3.551 1.884 + + lucid::vc(m3a) + ## effect component std.error z.ratio con + ## env!env.var 11.61 4.385 2.6 Positive + ## env:yorc!env.var 0.00063 0.000236 2.7 Positive + ## R!variance 3.551 0.2231 16 Positive + } + } } \keyword{datasets} diff --git a/man/piepho.barley.uniformity.Rd b/man/piepho.barley.uniformity.Rd index 71160ef5..651f757c 100644 --- a/man/piepho.barley.uniformity.Rd +++ b/man/piepho.barley.uniformity.Rd @@ -47,22 +47,23 @@ tick=TRUE, aspect=(36*3.73)/(30*1.90), main="piepho.barley.uniformity.csv") - libs(asreml,dplyr) - dat <- mutate(dat, x=factor(col), y=factor(row)) - dat <- arrange(dat, x, y) + if(require("asreml", quietly=TRUE)){ + libs(asreml,dplyr,lucid) + dat <- mutate(dat, x=factor(col), y=factor(row)) + dat <- arrange(dat, x, y) - # Piepho AR1xAR1 model (in random term, NOT residual) - m1 <- asreml(data=dat, - yield ~ 1, - random = ~ x + y + ar1(x):ar1(y), - residual = ~ units, - na.action=na.method(x="keep") ) - m1 <- update(m1) - # Match Piepho table 3, footnote 4: .9671, .9705 for col,row correlation - # Note these parameters are basically at the boundary of the parameter - # space. Questionable fit. - libs(lucid) - lucid::vc(m1) + # Piepho AR1xAR1 model (in random term, NOT residual) + m1 <- asreml(data=dat, + yield ~ 1, + random = ~ x + y + ar1(x):ar1(y), + residual = ~ units, + na.action=na.method(x="keep") ) + m1 <- update(m1) + # Match Piepho table 3, footnote 4: .9671, .9705 for col,row correlation + # Note these parameters are basically at the boundary of the parameter + # space. Questionable fit. + lucid::vc(m1) + } } } \keyword{datasets} diff --git a/man/piepho.cocksfoot.Rd b/man/piepho.cocksfoot.Rd index 3b11bd78..ead98294 100644 --- a/man/piepho.cocksfoot.Rd +++ b/man/piepho.cocksfoot.Rd @@ -45,33 +45,34 @@ libs(lattice) # Gaussian, not gamma distn densityplot(~date|year, data=dat, main="piepho.cocksfoot - heading date") - - libs(mumm) # The mumm package can reproduce Piepho's results - levelplot(date ~ year*gen, dat) - # note mp(random,fixed) - mod3 <- mumm(date ~ -1 + gen + (1|year) + mp(year, gen), dat) - - # Compare to Piepho table 3, "full maximum likelihood" - mod3$sigmas^2 # variances for year:gen, residual match - # year mp year:gen Residual - # 17.70287377 0.02944158 0.49024737 - - # mod3$par_fix # fixed genotypes match + if(require("mumm", quietly=TRUE)){ + libs(mumm) # The mumm package can reproduce Piepho's results + + levelplot(date ~ year*gen, dat) + # note mp(random,fixed) + mod3 <- mumm(date ~ -1 + gen + (1|year) + mp(year, gen), dat) + + # Compare to Piepho table 3, "full maximum likelihood" + mod3$sigmas^2 # variances for year:gen, residual match + # year mp year:gen Residual + # 17.70287377 0.02944158 0.49024737 + + # mod3$par_fix # fixed genotypes match - # mod3$sdreport # estim/stderr - # Estimate Std. Error - # nu 49.0393183 1.55038652 - # nu 42.0889493 1.67597832 - # nu 45.3411252 1.59818620 - # etc + # mod3$sdreport # estim/stderr + # Estimate Std. Error + # nu 49.0393183 1.55038652 + # nu 42.0889493 1.67597832 + # nu 45.3411252 1.59818620 + # etc - # mod3$par_rand # random year:gen match - # $`mp year:gen` - # 1990 1991 1992 1993 1994 1995 - # 0.10595661 -0.05298523 0.08228274 -0.09629696 -0.11045540 0.29637268 + # mod3$par_rand # random year:gen match + # $`mp year:gen` + # 1990 1991 1992 1993 1994 1995 + # 0.10595661 -0.05298523 0.08228274 -0.09629696 -0.11045540 0.29637268 + } } - } \keyword{datasets} diff --git a/man/shafi.tomato.uniformity.Rd b/man/shafi.tomato.uniformity.Rd index 6abd4a7a..a1739766 100644 --- a/man/shafi.tomato.uniformity.Rd +++ b/man/shafi.tomato.uniformity.Rd @@ -40,7 +40,7 @@ \dontrun{ library(agridat) data(shafi.tomato.uniformity) - shafi.tomato.uniformity <- dat + dat <- shafi.tomato.uniformity libs(desplot) desplot(dat, yield ~ col*row, diff --git a/man/sinclair.clover.Rd b/man/sinclair.clover.Rd index d31b8ce5..6d648a1e 100644 --- a/man/sinclair.clover.Rd +++ b/man/sinclair.clover.Rd @@ -97,7 +97,7 @@ title3d("sinclair.clover - yield","", xlab="Phosphorous/10", view3d(userMatrix=matrix(c(.7,.2,-.7,0, -.7,.2,-.6,0, 0,.9,.3,0, 0,0,0,1),ncol=4)) # snapshot3d(file, "png") -rgl.close() +close3d() } } diff --git a/man/smith.corn.uniformity.Rd b/man/smith.corn.uniformity.Rd index ff566228..bc7bbe8a 100644 --- a/man/smith.corn.uniformity.Rd +++ b/man/smith.corn.uniformity.Rd @@ -131,7 +131,7 @@ # black/gray dots very close to each other plot3d(dat$col, dat$row, dat$yield, col=dat$year, xlab="col",ylab="row",zlab="yield") - rgl.close() + close3d() } } diff --git a/man/smith.wheat.uniformity.Rd b/man/smith.wheat.uniformity.Rd index 734ba745..b87aa5e2 100644 --- a/man/smith.wheat.uniformity.Rd +++ b/man/smith.wheat.uniformity.Rd @@ -79,6 +79,7 @@ \dontrun{ library(agridat) + data(smith.wheat.uniformity) dat <- smith.wheat.uniformity libs(desplot) diff --git a/man/snedecor.asparagus.Rd b/man/snedecor.asparagus.Rd index 13444b68..f520dde7 100644 --- a/man/snedecor.asparagus.Rd +++ b/man/snedecor.asparagus.Rd @@ -61,63 +61,65 @@ # ---------- - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) - # Split-plot with asreml - m2 <- asreml(yield ~ trt + year + trt:year, data=dat, - random = ~ block + block:trt) - vc(m2) - ## effect component std.error z.ratio bound %ch - ## block 354.3 405 0.87 P 0.1 - ## block:trt 462.8 256.9 1.8 P 0 - ## units!R 404.7 82.6 4.9 P 0 + # Split-plot with asreml + m2 <- asreml(yield ~ trt + year + trt:year, data=dat, + random = ~ block + block:trt) + lucid::vc(m2) + ## effect component std.error z.ratio bound %ch + ## block 354.3 405 0.87 P 0.1 + ## block:trt 462.8 256.9 1.8 P 0 + ## units!R 404.7 82.6 4.9 P 0 + + ## # Antedependence with asreml. See O'Neill (2010). + dat <- dat[order(dat$block, dat$trt), ] + m3 <- asreml(yield ~ year * trt, data=dat, + random = ~ block, + residual = ~ block:trt:ante(year,1), + max=50) + m3 <- update(m3) + m3 <- update(m3) - ## # Antedependence with asreml. See O'Neill (2010). - dat <- dat[order(dat$block, dat$trt), ] - m3 <- asreml(yield ~ year * trt, data=dat, - random = ~ block, - residual = ~ block:trt:ante(year,1), - max=50) - m3 <- update(m3) - m3 <- update(m3) - - ## # Extract the covariance matrix for years and convert to correlation - ## covmat <- diag(4) - ## covmat[upper.tri(covmat,diag=TRUE)] <- m3$R.param$`block:trt:year`$year$initial - ## covmat[lower.tri(covmat)] <- t(covmat)[lower.tri(covmat)] - ## round(cov2cor(covmat),2) # correlation among the 4 years - ## # [,1] [,2] [,3] [,4] - ## # [1,] 1.00 0.45 0.39 0.31 - ## # [2,] 0.45 1.00 0.86 0.69 - ## # [3,] 0.39 0.86 1.00 0.80 - ## # [4,] 0.31 0.69 0.80 1.00 - - ## # We can also build the covariance Sigma by hand from the estimated - ## # variance components via: Sigma^-1 = U D^-1 U' - ## vv <- vc(m3) - ## print(vv) - ## ## effect component std.error z.ratio constr - ## ## block!block.var 86.56 156.9 0.55 pos - ## ## R!variance 1 NA NA fix - ## ## R!year.1930:1930 0.00233 0.00106 2.2 uncon - ## ## R!year.1931:1930 -0.7169 0.4528 -1.6 uncon - ## ## R!year.1931:1931 0.00116 0.00048 2.4 uncon - ## ## R!year.1932:1931 -1.139 0.1962 -5.8 uncon - ## ## R!year.1932:1932 0.00208 0.00085 2.4 uncon - ## ## R!year.1933:1932 -0.6782 0.1555 -4.4 uncon - ## ## R!year.1933:1933 0.00201 0.00083 2.4 uncon + ## # Extract the covariance matrix for years and convert to correlation + ## covmat <- diag(4) + ## covmat[upper.tri(covmat,diag=TRUE)] <- m3$R.param$`block:trt:year`$year$initial + ## covmat[lower.tri(covmat)] <- t(covmat)[lower.tri(covmat)] + ## round(cov2cor(covmat),2) # correlation among the 4 years + ## # [,1] [,2] [,3] [,4] + ## # [1,] 1.00 0.45 0.39 0.31 + ## # [2,] 0.45 1.00 0.86 0.69 + ## # [3,] 0.39 0.86 1.00 0.80 + ## # [4,] 0.31 0.69 0.80 1.00 + + ## # We can also build the covariance Sigma by hand from the estimated + ## # variance components via: Sigma^-1 = U D^-1 U' + ## vv <- vc(m3) + ## print(vv) + ## ## effect component std.error z.ratio constr + ## ## block!block.var 86.56 156.9 0.55 pos + ## ## R!variance 1 NA NA fix + ## ## R!year.1930:1930 0.00233 0.00106 2.2 uncon + ## ## R!year.1931:1930 -0.7169 0.4528 -1.6 uncon + ## ## R!year.1931:1931 0.00116 0.00048 2.4 uncon + ## ## R!year.1932:1931 -1.139 0.1962 -5.8 uncon + ## ## R!year.1932:1932 0.00208 0.00085 2.4 uncon + ## ## R!year.1933:1932 -0.6782 0.1555 -4.4 uncon + ## ## R!year.1933:1933 0.00201 0.00083 2.4 uncon + + ## U <- diag(4) + ## U[1,2] <- vv[4,2] ; U[2,3] <- vv[6,2] ; U[3,4] <- vv[8,2] + ## Dinv <- diag(c(vv[3,2], vv[5,2], vv[7,2], vv[9,2])) + ## # solve(U %*% Dinv %*% t(U)) # same as 'covmat' above + ## solve(crossprod(t(U), tcrossprod(Dinv, U)) ) + ## ## [,1] [,2] [,3] [,4] + ## ## [1,] 428.4310 307.1478 349.8152 237.2453 + ## ## [2,] 307.1478 1083.9717 1234.5516 837.2751 + ## ## [3,] 349.8152 1234.5516 1886.5150 1279.4378 + ## ## [4,] 237.2453 837.2751 1279.4378 1364.8446 + } - ## U <- diag(4) - ## U[1,2] <- vv[4,2] ; U[2,3] <- vv[6,2] ; U[3,4] <- vv[8,2] - ## Dinv <- diag(c(vv[3,2], vv[5,2], vv[7,2], vv[9,2])) - ## # solve(U %*% Dinv %*% t(U)) # same as 'covmat' above - ## solve(crossprod(t(U), tcrossprod(Dinv, U)) ) - ## ## [,1] [,2] [,3] [,4] - ## ## [1,] 428.4310 307.1478 349.8152 237.2453 - ## ## [2,] 307.1478 1083.9717 1234.5516 837.2751 - ## ## [3,] 349.8152 1234.5516 1886.5150 1279.4378 - ## ## [4,] 237.2453 837.2751 1279.4378 1364.8446 - } } \keyword{datasets} diff --git a/man/steptoe.morex.pheno.Rd b/man/steptoe.morex.pheno.Rd index 2fea4dee..7d20aaeb 100644 --- a/man/steptoe.morex.pheno.Rd +++ b/man/steptoe.morex.pheno.Rd @@ -159,31 +159,35 @@ libs(qtl) data(steptoe.morex.geno) datg <- steptoe.morex.geno - plot.map(datg, main="steptoe.morex.pheno") # or just use plot() - + qtl::plot.map(datg, main="steptoe.morex.pheno") + qtl::plotMissing(datg) + # This is a very rudimentary example. + # The 'wgaim' function works interactively, but fails during + # devtools::check(). + if(0 & require("asreml", quietly=TRUE)){ + libs(asreml) + + # Fit a simple multi-environment mixed model + m1 <- asreml(yield ~ env, data=dat, random=~gen) + + libs(wgaim) + wgaim::linkMap(datg) + # Create an interval object for wgaim + dati <- wgaim::cross2int(datg, id="gen") - # Fit a simple multi-environment mixed model - libs(asreml) - m1 <- asreml(yield ~ env, data=dat, random=~gen) - - libs(wgaim) - qtl::plotMissing(datg) - wgaim::linkMap(datg) - # Create an interval object for wgaim - dati <- wgaim::cross2int(datg, id="gen") - - # Whole genome qtl - q1 <- wgaim::wgaim(m1, intervalObj=dati, merge.by="gen", na.action=na.method(x="include")) - #wgaim::linkMap(q1, dati) # Visualize - wgaim::outStat(q1, dati) # outlier statistic - summary(q1, dati) # Table of important intervals - # Chrom Left Marker dist(cM) Right Marker dist(cM) Size Pvalue % Var - # 3 ABG399 52.6 BCD828 56.1 0.254 0.000 45.0 - # 5 MWG912 148 ABG387A 151.2 0.092 0.001 5.9 - # 6 ABC169B 64.8 CDO497 67.5 -0.089 0.001 5.6 - + # Whole genome qtl + q1 <- wgaim::wgaim(m1, intervalObj=dati, + merge.by="gen", na.action=na.method(x="include")) + #wgaim::linkMap(q1, dati) # Visualize + wgaim::outStat(q1, dati) # outlier statistic + summary(q1, dati) # Table of important intervals + # Chrom Left Marker dist(cM) Right Marker dist(cM) Size Pvalue % Var + # 3 ABG399 52.6 BCD828 56.1 0.254 0.000 45.0 + # 5 MWG912 148 ABG387A 151.2 0.092 0.001 5.9 + # 6 ABC169B 64.8 CDO497 67.5 -0.089 0.001 5.6 + } } } diff --git a/man/stroup.nin.Rd b/man/stroup.nin.Rd index 93557427..704ddd9d 100644 --- a/man/stroup.nin.Rd +++ b/man/stroup.nin.Rd @@ -177,36 +177,36 @@ # ----- - # asreml - libs(asreml,lucid) - - # RCB analysis - as1 <- asreml(yield ~ gen, random = ~ rep, data=dat, - na.action=na.method(x="omit")) - preds$asreml1 <- predict(as1, data=dat, classify="gen")$pvals$predicted.value - - # Two-dimensional AR1xAR1 spatial model - dat <- transform(dat, xf=factor(col), yf=factor(row)) - dat <- dat[order(dat$xf, dat$yf),] - as2 <- asreml(yield~gen, data=dat, - residual = ~ar1(xf):ar1(yf), - na.action=na.method(x="omit")) - preds$asreml2 <- predict(as2, data=dat, classify="gen")$pvals$predicted.value - - lucid::vc(as2) - ## effect component std.error z.ratio constr - ## R!variance 48.7 7.155 6.8 pos - ## R!xf.cor 0.6555 0.05638 12 unc - ## R!yf.cor 0.4375 0.0806 5.4 unc + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) + + # RCB analysis + as1 <- asreml(yield ~ gen, random = ~ rep, data=dat, + na.action=na.method(x="omit")) + preds$asreml1 <- predict(as1, data=dat, classify="gen")$pvals$predicted.value + + # Two-dimensional AR1xAR1 spatial model + dat <- transform(dat, xf=factor(col), yf=factor(row)) + dat <- dat[order(dat$xf, dat$yf),] + as2 <- asreml(yield~gen, data=dat, + residual = ~ar1(xf):ar1(yf), + na.action=na.method(x="omit")) + preds$asreml2 <- predict(as2, data=dat, classify="gen")$pvals$predicted.value + + lucid::vc(as2) + ## effect component std.error z.ratio constr + ## R!variance 48.7 7.155 6.8 pos + ## R!xf.cor 0.6555 0.05638 12 unc + ## R!yf.cor 0.4375 0.0806 5.4 unc # Compare the estimates from the two asreml models. # We see that Buckskin has correctly been shifted upward by the spatial model - plot(preds$as1, preds$as2, xlim=c(13,37), ylim=c(13,37), - xlab="RCB", ylab="AR1xAR1", type='n') - title("stroup.nin: Comparison of predicted values") - text(preds$asreml1, preds$asreml2, preds$gen, cex=0.5) - abline(0,1) - + plot(preds$as1, preds$as2, xlim=c(13,37), ylim=c(13,37), + xlab="RCB", ylab="AR1xAR1", type='n') + title("stroup.nin: Comparison of predicted values") + text(preds$asreml1, preds$asreml2, preds$gen, cex=0.5) + abline(0,1) + } # ----- # sommer diff --git a/man/stroup.splitplot.Rd b/man/stroup.splitplot.Rd index 973cfe6f..3fd29504 100644 --- a/man/stroup.splitplot.Rd +++ b/man/stroup.splitplot.Rd @@ -56,38 +56,38 @@ # m0 <- lme(y ~ -1 + a + b + a:b, data=dat, random = ~ 1|rep/a) # ----- ASREML model --- - libs(asreml) - m1 <- asreml(y~ -1 + a + b + a:b, random=~ rep + a:rep, data=dat) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) + m1 <- asreml(y~ -1 + a + b + a:b, random=~ rep + a:rep, data=dat) - libs(lucid) - # vc(m1) # Variance components match Stroup p. 41 - ## effect component std.error z.ratio bound - ## rep 62.42 56.41 1.1 P - ## a:rep 15.39 11.8 1.3 P - ## units(R) 9.364 4.415 2.1 P + # vc(m1) # Variance components match Stroup p. 41 + ## effect component std.error z.ratio bound + ## rep 62.42 56.41 1.1 P + ## a:rep 15.39 11.8 1.3 P + ## units(R) 9.364 4.415 2.1 P - # Narrow space predictions - predict(m1, data=dat, classify="a", average=list(rep=NULL)) - # a Predicted Std Err Status - # a1 32.88 1.082 Estimable - # a2 34.12 1.082 Estimable - # a3 25.75 1.082 Estimable - - # Intermediate space predictions - predict(m1, data=dat, classify="a", ignore="a:rep", - average=list(rep=NULL)) - # a Predicted Std Err Status - # a1 32.88 2.24 Estimable - # a2 34.12 2.24 Estimable - # a3 25.75 2.24 Estimable - - # Broad space predictions - predict(m1, data=dat, classify="a") - # a Predicted Std Err Status - # a1 32.88 4.54 Estimable - # a2 34.12 4.54 Estimable - # a3 25.75 4.54 Estimable - + # Narrow space predictions + predict(m1, data=dat, classify="a", average=list(rep=NULL)) + # a Predicted Std Err Status + # a1 32.88 1.082 Estimable + # a2 34.12 1.082 Estimable + # a3 25.75 1.082 Estimable + + # Intermediate space predictions + predict(m1, data=dat, classify="a", ignore="a:rep", + average=list(rep=NULL)) + # a Predicted Std Err Status + # a1 32.88 2.24 Estimable + # a2 34.12 2.24 Estimable + # a3 25.75 2.24 Estimable + + # Broad space predictions + predict(m1, data=dat, classify="a") + # a Predicted Std Err Status + # a1 32.88 4.54 Estimable + # a2 34.12 4.54 Estimable + # a3 25.75 4.54 Estimable + } # ----- MCMCglmm model ----- # Use the point estimates from REML with a prior distribution @@ -115,7 +115,7 @@ # the linear contrast for the the narrow-space predictions # (also called adjusted mean) for the a1:a2 interaction. # a1n a1i a1b a1a2n a1a2ib - cm <- matrix(c(1, 1, 1, 1, 1, # a1 + cm <- matrix(c(1, 1, 1, 1, 1, # a1 0, 0, 0, -1, -1, # a2 0, 0, 0, 0, 0, # a3 1/2, 1/2, 1/2, 0, 0, # b2 @@ -159,13 +159,15 @@ # REML est 32.88 32.88 32.88 -1.25 -1.25 # REML stderr 1.08 2.24 4.54 1.53 3.17 # MCMC mode 32.95 32.38 31.96 -1.07 -1.17 - # MCMC stderr 1.23 2.64 5.93 1.72 3.73 - plot(post2) - + # MCMC stderr 1.23 2.64 5.93 1.72 3.73 + # plot(post2) + post22 <- lattice::make.groups( Narrow=post2[,1], Intermediate=post2[,2], Broad=post2[,3]) print(densityplot(~data|which, data=post22, groups=which, cex=.25, lty=1, layout=c(1,3), + main="stroup.splitplot", xlab="MCMC model value of predictable function for A1")) + } } diff --git a/man/tesfaye.millet.Rd b/man/tesfaye.millet.Rd index aae797ef..3b9f9806 100644 --- a/man/tesfaye.millet.Rd +++ b/man/tesfaye.millet.Rd @@ -26,7 +26,7 @@ Experiments conducted at Bako and Assosa research centers in Ethiopia. The data has: - 4 year, 2 sites = 7 environments, + 4 years, 2 sites = 7 environments, 2-3 reps per trial, 47 genotypes. @@ -49,30 +49,33 @@ \examples{ \dontrun{ -library(agridat) -data(tesfaye.millet) -dat <- tesfaye.millet + library(agridat) + data(tesfaye.millet) + dat <- tesfaye.millet -libs(desplot) -desplot(dat, yield ~ col*row|env, out1=rep, main="tesfaye.millet") + dat <- transform(dat, year=factor(year), site=factor(site)) + libs(dplyr,asreml,lucid) + dat <- mutate(dat, + env=factor(paste0(site,year)), + gen=factor(gen), + rep=factor(rep), + xfac=factor(col), yfac=factor(row)) + libs(desplot) + desplot(dat, yield~col*row|env, main="tesfaye.millet") + dat <- arrange(dat, env, xfac, yfac) -libs(dplyr,asreml,lucid) -dat <- mutate(dat, - env=factor(paste0(site,year)), - gen=factor(gen), - rep=factor(rep), - xfac=factor(col), yfac=factor(row)) -dat <- arrange(dat, env, xfac, yfac) - -# Fixed environment -# Random row/col within environment, Factor Analytic GxE -# AR1xAR1 spatial residuals within each environment -m1 <- asreml(yield ~ 1 + env, - data=dat, - random = ~ at(env):xfac + at(env):yfac + gen:fa(env), - residual = ~ dsum( ~ ar1(xfac):ar1(yfac)|env) ) -m1 <- update(m1) -lucid::vc(m1) + # Fixed environment + # Random row/col within environment, Factor Analytic GxE + # AR1xAR1 spatial residuals within each environment + if(require("asreml", quietly=TRUE)){ + libs(asreml) + m1 <- asreml(yield ~ 1 + env, + data=dat, + random = ~ at(env):xfac + at(env):yfac + gen:fa(env), + residual = ~ dsum( ~ ar1(xfac):ar1(yfac)|env) ) + m1 <- update(m1) + lucid::vc(m1) + } } } diff --git a/man/thompson.cornsoy.Rd b/man/thompson.cornsoy.Rd index 02912ecc..59e41344 100644 --- a/man/thompson.cornsoy.Rd +++ b/man/thompson.cornsoy.Rd @@ -122,7 +122,7 @@ xyplot(corn+soy~year|state, dat, m2 <- aov(corn ~ ns(rain0, 3) + ns(rain7, 3) + ns(temp8, 3) + ns(year,3), dat.ia) op <- par(mfrow=c(2,2)) - termplot(m2, se=TRUE, rug=TRUE, partial=TRUE) + termplot(m2, se=TRUE, rug=TRUE, partial=TRUE, main="thompson.cornsoy") par(op) # do NOT use gam package @@ -130,7 +130,7 @@ xyplot(corn+soy~year|state, dat, m1 <- gam(corn ~ s(year, k=5) + s(rain0, k=5) + s(rain7, k=5) + s(temp8, k=5), data=dat.ia) op <- par(mfrow=c(2,2)) - plot.gam(m1, residuals=TRUE, se=TRUE, cex=2) + plot.gam(m1, residuals=TRUE, se=TRUE, cex=2, main="thompson.cornsoy") par(op) } diff --git a/man/usgs.herbicides.Rd b/man/usgs.herbicides.Rd index 2e586587..b623e72c 100644 --- a/man/usgs.herbicides.Rd +++ b/man/usgs.herbicides.Rd @@ -76,7 +76,8 @@ Sample types: CR = concurrent replicate sample, FB = field blank, LD = laborator # median/mean with(dat, cenmle(y, ycen, dist="lognormal")) # boxplot - with(dat, cenboxplot(obs=y, cen=ycen, log=FALSE)) + with(dat, cenboxplot(obs=y, cen=ycen, log=FALSE, + main="usgs.herbicides")) # with(dat, boxplot(y)) pp <- with(dat, ros(obs=y, censored=ycen, forwardT="log")) # default lognormal plot(pp) diff --git a/man/verbyla.lupin.Rd b/man/verbyla.lupin.Rd index 8727807c..6a22fcb2 100644 --- a/man/verbyla.lupin.Rd +++ b/man/verbyla.lupin.Rd @@ -117,88 +117,90 @@ # ---------- - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) - # We try to reproduce the analysis of Verbyla 1999. - # May not be exactly the same, but is pretty close. + # We try to reproduce the analysis of Verbyla 1999. + # May not be exactly the same, but is pretty close. - # Check nlevels for size of random-coefficient structures - # length(with(dat, table(gen))) # 9 varieties for RC1 - # length(with(dat, table(gen,site))) # 99 site:gen combinations for RC2 - - # Make row and col into factors - dat <- transform(dat, colf=factor(col), rowf=factor(row)) - # sort for asreml - dat <- dat[order(dat$site, dat$rowf, dat$colf),] + # Check nlevels for size of random-coefficient structures + # length(with(dat, table(gen))) # 9 varieties for RC1 + # length(with(dat, table(gen,site))) # 99 site:gen combinations for RC2 - # Make site names more useful for plots - # dat <- transform(dat, site=factor(paste0(year,".",substring(loc,1,4)))) + # Make row and col into factors + dat <- transform(dat, colf=factor(col), rowf=factor(row)) + # sort for asreml + dat <- dat[order(dat$site, dat$rowf, dat$colf),] + + # Make site names more useful for plots + # dat <- transform(dat, site=factor(paste0(year,".",substring(loc,1,4)))) + + # Initial model from top of Verbyla table 9. + m0 <- asreml(yield ~ 1 + + site + + linrate + + site:linrate, + data = dat, + random = ~ spl(rate) + + dev(rate) + + site:spl(rate) + + site:dev(rate) + + str(~gen+gen:linrate, ~us(2):id(9)) # RC1 + + gen:spl(rate) + + gen:dev(rate) + + str(~site:gen+site:gen:linrate, ~us(2):id(99)) # RC2 + + site:gen:spl(rate) + + site:gen:dev(rate), + residual = ~ dsum( ~ ar1(rowf):ar1(colf)|site) # Spatial AR1 x AR1 + ) + m0 <- update(m0) + m0 <- update(m0) + m0 <- update(m0) + m0 <- update(m0) + m0 <- update(m0) + + # Variograms match Verbyla 1999 figure 7 (scale slightly different) + plot(varioGram(m0), xlim=c(1:19), zlim=c(0,2), + main="verbyla.lupin - variogram by site") - # Initial model from top of Verbyla table 9. - m0 <- asreml(yield ~ 1 - + site - + linrate - + site:linrate, - data = dat, - random = ~ spl(rate) - + dev(rate) - + site:spl(rate) - + site:dev(rate) - + str(~gen+gen:linrate, ~us(2):id(9)) # RC1 - + gen:spl(rate) - + gen:dev(rate) - + str(~site:gen+site:gen:linrate, ~us(2):id(99)) # RC2 - + site:gen:spl(rate) - + site:gen:dev(rate), - residual = ~ dsum( ~ ar1(rowf):ar1(colf)|site) # Spatial AR1 x AR1 - ) - m0 <- update(m0) - m0 <- update(m0) - m0 <- update(m0) - m0 <- update(m0) - m0 <- update(m0) - - # Variograms match Verbyla 1999 figure 7 (scale slightly different) - plot(varioGram(m0), xlim=c(1:19), zlim=c(0,2), - main="verbyla.lupin - variogram by site") - # Sequence of models in Verbyla 1999 table 10 - m1 <- update(m0, fixed= ~ . - + at(site, c(2,5,6,8,9,10)):lincol - + at(site, c(3,5,7,8)):linrow - + at(site, c(2,3,5,7,8,9,11)):serp + m1 <- update(m0, fixed= ~ . + + at(site, c(2,5,6,8,9,10)):lincol + + at(site, c(3,5,7,8)):linrow + + at(site, c(2,3,5,7,8,9,11)):serp , random = ~ . - + at(site, c(3,6,7,9)):rowf - + at(site, c(1,2,3,9,10)):colf - + at(site, c(5,7,8,10)):units) - m1 <- update(m1) - - m2 <- update(m1, - random = ~ . - - site:gen:spl(rate) - site:gen:dev(rate)) - - m3 <- update(m2, - random = ~ . - - site:dev(rate) - gen:dev(rate)) - - m4 <- update(m3, - random = ~ . - - dev(rate)) - - m5 <- update(m4, - random = ~ . - - at(site, c(5,7,8,10)):units + at(site, c(5,7,8)):units) - - # Variance components are a pretty good match to Verbyla 1997, table 15 - libs(lucid) - vc(m5) - .001004/sqrt(.005446*.0003662) # .711 correlation for RC1 - .00175/sqrt(.01881*.000167) # .987 correlation for RC2 - - # Matches Verbyla 1999 figure 5 - plot(varioGram(m5), - main="verbyla.lupin - final model variograms", - xlim=c(1:19), zlim=c(0,1.5)) + + at(site, c(3,6,7,9)):rowf + + at(site, c(1,2,3,9,10)):colf + + at(site, c(5,7,8,10)):units) + m1 <- update(m1) + + m2 <- update(m1, + random = ~ . + - site:gen:spl(rate) - site:gen:dev(rate)) + + m3 <- update(m2, + random = ~ . + - site:dev(rate) - gen:dev(rate)) + + m4 <- update(m3, + random = ~ . + - dev(rate)) + + m5 <- update(m4, + random = ~ . + - at(site, c(5,7,8,10)):units + at(site, c(5,7,8)):units) + + # Variance components are a pretty good match to Verbyla 1997, table 15 + libs(lucid) + vc(m5) + .001004/sqrt(.005446*.0003662) # .711 correlation for RC1 + .00175/sqrt(.01881*.000167) # .987 correlation for RC2 + + # Matches Verbyla 1999 figure 5 + plot(varioGram(m5), + main="verbyla.lupin - final model variograms", + xlim=c(1:19), zlim=c(0,1.5)) + } } } diff --git a/man/vold.longterm.Rd b/man/vold.longterm.Rd index cb903531..6bdfc0bf 100644 --- a/man/vold.longterm.Rd +++ b/man/vold.longterm.Rd @@ -78,7 +78,8 @@ # Separate fit for each year. Overfitting with 3x19=57 params. libs(nlme) m2.lis <- nlsList(yield ~ SSasymp(nitro,max,int,lograte) | year, data=dat) - plot(intervals(m2.lis),layout = c(3,1)) # lograte might be same for each year + plot(intervals(m2.lis),layout = c(3,1), + main="vold.longterm") # lograte might be same for each year # Fixed overall asymptotic model, plus random deviations for each year diff --git a/man/vsn.lupin3.Rd b/man/vsn.lupin3.Rd index c07fb036..96933258 100644 --- a/man/vsn.lupin3.Rd +++ b/man/vsn.lupin3.Rd @@ -55,92 +55,92 @@ desplot(dat, check~ col*row|site, main="vsn.lupin3: check plot placement") - # asreml - libs(asreml,lucid) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) - # Single-site analyses suggested random row term for site 3, - # random column terms for all sites, - # AR1 was unnecessary for the col dimension of site 3 - dat <- transform(dat, colf=factor(col), rowf=factor(row)) - dat <- dat[order(dat$site, dat$colf, dat$rowf),] # Sort for asreml - m1 <- asreml(yield ~ site + check:site, data=dat, - random = ~ at(site):colf + at(site,3):rowf + test, - residual = ~ dsum( ~ ar1(colf):ar1(rowf) + - id(colf):ar1(rowf) | site, - levels=list(1:2, 3) - ) ) - m1$loglik - ## [1] -314.2616 + # Single-site analyses suggested random row term for site 3, + # random column terms for all sites, + # AR1 was unnecessary for the col dimension of site 3 + dat <- transform(dat, colf=factor(col), rowf=factor(row)) + dat <- dat[order(dat$site, dat$colf, dat$rowf),] # Sort for asreml + m1 <- asreml(yield ~ site + check:site, data=dat, + random = ~ at(site):colf + at(site,3):rowf + test, + residual = ~ dsum( ~ ar1(colf):ar1(rowf) + + id(colf):ar1(rowf) | site, + levels=list(1:2, 3) + ) ) + m1$loglik + ## [1] -314.2616 + + lucid::vc(m1) + ## effect component std.error z.ratio constr + ## at(site, S1):colf!colf.var 0.6228 0.4284 1.5 pos + ## at(site, S2):colf!colf.var 0.159 0.1139 1.4 pos + ## at(site, S3):colf!colf.var 0.04832 0.02618 1.8 pos + ## at(site, S3):rowf!rowf.var 0.0235 0.008483 2.8 pos + ## test!test.var 0.1031 0.01468 7 pos + ## site_S1!variance 2.771 0.314 8.8 pos + ## site_S1!colf.cor 0.1959 0.05375 3.6 uncon + ## site_S1!rowf.cor 0.6503 0.03873 17 uncon + ## site_S2!variance 0.9926 0.1079 9.2 pos + ## site_S2!colf.cor 0.2868 0.05246 5.5 uncon + ## site_S2!rowf.cor 0.5744 0.0421 14 uncon + ## site_S3!variance 0.1205 0.01875 6.4 pos + ## site_S3!rowf.cor 0.6394 0.06323 10 uncon + + # Add site:test + m2 <- update(m1, random=~. + site:test) + m2$loglik + ## [1] -310.8794 + + # CORUH structure on the site component of site:test + m3 <- asreml(yield ~ site + check:site, data=dat, + random = ~ at(site):colf + at(site,3):rowf + corh(site):test, + residual = ~ dsum( ~ ar1(colf):ar1(rowf) + + id(colf):ar1(rowf) | site, + levels=list(1:2, 3) )) + m3$loglik + ## [1] -288.4837 + + # Unstructured genetic variance matrix + m4 <- asreml(yield ~ site + check:site, data=dat, + random = ~ at(site):colf + at(site,3):rowf + us(site):test, + residual = ~ dsum( ~ ar1(colf):ar1(rowf) + + id(colf):ar1(rowf) | site, + levels=list(1:2, 3) )) + m4$loglik + ## [1] -286.8239 - vc(m1) - ## effect component std.error z.ratio constr - ## at(site, S1):colf!colf.var 0.6228 0.4284 1.5 pos - ## at(site, S2):colf!colf.var 0.159 0.1139 1.4 pos - ## at(site, S3):colf!colf.var 0.04832 0.02618 1.8 pos - ## at(site, S3):rowf!rowf.var 0.0235 0.008483 2.8 pos - ## test!test.var 0.1031 0.01468 7 pos - ## site_S1!variance 2.771 0.314 8.8 pos - ## site_S1!colf.cor 0.1959 0.05375 3.6 uncon - ## site_S1!rowf.cor 0.6503 0.03873 17 uncon - ## site_S2!variance 0.9926 0.1079 9.2 pos - ## site_S2!colf.cor 0.2868 0.05246 5.5 uncon - ## site_S2!rowf.cor 0.5744 0.0421 14 uncon - ## site_S3!variance 0.1205 0.01875 6.4 pos - ## site_S3!rowf.cor 0.6394 0.06323 10 uncon - - # Add site:test - m2 <- update(m1, random=~. + site:test) - m2$loglik - ## [1] -310.8794 - - # CORUH structure on the site component of site:test - m3 <- asreml(yield ~ site + check:site, data=dat, - random = ~ at(site):colf + at(site,3):rowf + corh(site):test, - residual = ~ dsum( ~ ar1(colf):ar1(rowf) + - id(colf):ar1(rowf) | site, - levels=list(1:2, 3) )) - m3$loglik - ## [1] -288.4837 - - # Unstructured genetic variance matrix - m4 <- asreml(yield ~ site + check:site, data=dat, - random = ~ at(site):colf + at(site,3):rowf + us(site):test, - residual = ~ dsum( ~ ar1(colf):ar1(rowf) + - id(colf):ar1(rowf) | site, - levels=list(1:2, 3) )) - m4$loglik - ## [1] -286.8239 - - # Note that a 3x3 unstructured matrix can be written LL'+Psi with 1 factor L - # Explicitly fit the factor analytic model - m5 <- asreml(yield ~ site + check:site, data=dat, - random = ~ at(site):colf + at(site,3):rowf - + fa(site,1, init=c(.7,.1,.1,.5,.3,.2)):test, - residual = ~ dsum( ~ ar1(colf):ar1(rowf) + - id(colf):ar1(rowf) | site, - levels=list(1:2, 3) )) - m5$loglik # Same as m4 - ## [1] -286.8484 - - # Model 4, Unstructured (symmetric) genetic variance matrix - un <- diag(3) - un[upper.tri(un,TRUE)] <- m4$vparameters[5:10] - round(un+t(un)-diag(diag(un)),3) - ## [,1] [,2] [,3] - ## [1,] 0.992 0.158 0.132 - ## [2,] 0.158 0.073 0.078 - ## [3,] 0.132 0.078 0.122 - - # Model 5, FA matrix = LL'+Psi. Not quite the same as unstructured, - # since the FA model fixes site 2 variance at 0. - psi <- diag(m5$vparameters[5:7]) - lam <- matrix(m5$vparameters[8:10], ncol=1) - round(tcrossprod(lam,lam)+psi,3) - ## [,1] [,2] [,3] - ## [1,] 0.991 0.156 0.133 - ## [2,] 0.156 0.092 0.078 - ## [3,] 0.133 0.078 0.122 - + # Note that a 3x3 unstructured matrix can be written LL'+Psi with 1 factor L + # Explicitly fit the factor analytic model + m5 <- asreml(yield ~ site + check:site, data=dat, + random = ~ at(site):colf + at(site,3):rowf + + fa(site,1, init=c(.7,.1,.1,.5,.3,.2)):test, + residual = ~ dsum( ~ ar1(colf):ar1(rowf) + + id(colf):ar1(rowf) | site, + levels=list(1:2, 3) )) + m5$loglik # Same as m4 + ## [1] -286.8484 + + # Model 4, Unstructured (symmetric) genetic variance matrix + un <- diag(3) + un[upper.tri(un,TRUE)] <- m4$vparameters[5:10] + round(un+t(un)-diag(diag(un)),3) + ## [,1] [,2] [,3] + ## [1,] 0.992 0.158 0.132 + ## [2,] 0.158 0.073 0.078 + ## [3,] 0.132 0.078 0.122 + + # Model 5, FA matrix = LL'+Psi. Not quite the same as unstructured, + # since the FA model fixes site 2 variance at 0. + psi <- diag(m5$vparameters[5:7]) + lam <- matrix(m5$vparameters[8:10], ncol=1) + round(tcrossprod(lam,lam)+psi,3) + ## [,1] [,2] [,3] + ## [1,] 0.991 0.156 0.133 + ## [2,] 0.156 0.092 0.078 + ## [3,] 0.133 0.078 0.122 + } } } diff --git a/man/wedderburn.barley.Rd b/man/wedderburn.barley.Rd index 071c6566..aa47668b 100644 --- a/man/wedderburn.barley.Rd +++ b/man/wedderburn.barley.Rd @@ -86,11 +86,12 @@ libs(gnm) m3 <- glm(y ~ gen + site, data=dat, family="wedderburn") summary(m3) # Similar to McCullagh fig 9.2 -plot(m3, which=1) + plot(m3, which=1) + title("wedderburn.barley") # Compare data and model dat$pwed <- predict(m3, type="response") -dotplot(gen~pwed+y|site, dat) +dotplot(gen~pwed+y|site, dat, main="wedderburn.barley") } } diff --git a/man/weiss.incblock.Rd b/man/weiss.incblock.Rd index daca69b1..4621b61e 100644 --- a/man/weiss.incblock.Rd +++ b/man/weiss.incblock.Rd @@ -60,18 +60,19 @@ main="weiss.incblock") - # asreml - # Standard inc block analysis used by Weiss and Cox - libs(asreml) - m1 <- asreml(yield ~ gen + block , data=dat) - predict(m1, data=dat, classify="gen")$pvals - - ## gen pred.value std.error est.stat - ## G01 24.59 0.8312 Estimable - ## G02 26.92 0.8312 Estimable - ## G03 32.62 0.8312 Estimable - ## G04 26.97 0.8312 Estimable - ## G05 26.02 0.8312 Estimable + if(require("asreml", quietly=TRUE)){ + # Standard inc block analysis used by Weiss and Cox + libs(asreml) + m1 <- asreml(yield ~ gen + block , data=dat) + predict(m1, data=dat, classify="gen")$pvals + + ## gen pred.value std.error est.stat + ## G01 24.59 0.8312 Estimable + ## G02 26.92 0.8312 Estimable + ## G03 32.62 0.8312 Estimable + ## G04 26.97 0.8312 Estimable + ## G05 26.02 0.8312 Estimable + } } } diff --git a/man/weiss.lattice.Rd b/man/weiss.lattice.Rd index 3b58a6bb..69adc943 100644 --- a/man/weiss.lattice.Rd +++ b/man/weiss.lattice.Rd @@ -65,20 +65,21 @@ rep (instead of across reps). # ---------- - libs(asreml) - m2 <- asreml(yield ~ rep + rep:xf + rep:yf + gen, data=dat) - # Weiss table 6 means - wald(m2) - predict(m2, data=dat, classify="gen")$pvals + if(require("asreml", quietly=TRUE)){ + libs(asreml) + m2 <- asreml(yield ~ rep + rep:xf + rep:yf + gen, data=dat) + # Weiss table 6 means + wald(m2) + predict(m2, data=dat, classify="gen")$pvals + ## gen pred.value std.error est.stat + ## G01 27.74 1.461 Estimable + ## G02 24.95 1.461 Estimable + ## G03 24.38 1.461 Estimable + ## G04 28.05 1.461 Estimable + ## G05 19.6 1.461 Estimable + ## G06 23.79 1.461 Estimable + } - ## gen pred.value std.error est.stat - ## G01 27.74 1.461 Estimable - ## G02 24.95 1.461 Estimable - ## G03 24.38 1.461 Estimable - ## G04 28.05 1.461 Estimable - ## G05 19.6 1.461 Estimable - ## G06 23.79 1.461 Estimable - } } \keyword{datasets} diff --git a/man/wiebe.wheat.uniformity.Rd b/man/wiebe.wheat.uniformity.Rd index 95720e68..65c574eb 100644 --- a/man/wiebe.wheat.uniformity.Rd +++ b/man/wiebe.wheat.uniformity.Rd @@ -132,6 +132,7 @@ # and found 21.5 g / inch of space between rows. We used all the # data without correcting for fertility and obtained 33.1 g / inch. xyplot(yield ~ area, dat, type=c('p','r'), + main="wiebe.wheat.uniformity", xlab="Average area per row", ylab="Yield") coef(lm(yield ~ area, dat))[2] # 33.1 diff --git a/man/yates.oats.Rd b/man/yates.oats.Rd index aa33049c..cd5c3aeb 100644 --- a/man/yates.oats.Rd +++ b/man/yates.oats.Rd @@ -51,8 +51,7 @@ \references{ Yates, Frank (1935) Complex experiments, - \emph{Journal of the Royal Statistical Society Suppl} - 2, 181-247. Figure 2. + Journal of the Royal Statistical Society Supplement 2, 181-247. Figure 2. https://doi.org/10.2307/2983638 } @@ -63,7 +62,6 @@ data(yates.oats) dat <- yates.oats - ## # Means match Rothamsted report p. 144 ## libs(dplyr) ## dat %>% group_by(nitro,gen) %>% @@ -105,7 +103,7 @@ # ---------- - # Marginal predictions from emmeans package and asreml::predict + # Marginal predictions # --- nlme --- libs(nlme) @@ -119,44 +117,17 @@ dat2$nitrof <- factor(dat2$nitro) # --- asreml --- - libs(asreml) - m5a <- asreml(yield ~ nitrof + gen, - random = ~ block + block:gen, data=dat2) - libs(lucid) - vc(m5l) - vc(m5a) + if(require("asreml", quietly=TRUE)){ + libs(asreml,lucid) + m5a <- asreml(yield ~ nitrof + gen, + random = ~ block + block:gen, data=dat2) + lucid::vc(m5l) + lucid::vc(m5a) emmeans::emmeans(m5l, "gen") predict(m5a, data=dat2, classify="gen")$pvals - - # ---------- - - if(0){ - - # Demonstrate use of regress package, compare to lme - - libs(regress) - m6 <- regress(yield ~ nitrof + gen, ~block + I(block:gen), identity=TRUE, - verbose=1, data=dat) - summary(m6) - ## Variance Coefficients: - ## Estimate Std. Error - ## block 214.468 168.794 - ## I(block:gen) 109.700 67.741 - ## In 162.558 32.189 - - # ordinal causes clash with VarCorr - if(is.element("package:ordinal", search())) detach(package:ordinal) - - m7 <- lme(yield ~ nitrof + gen, random = ~ 1|block/gen, data=dat) - lme4::VarCorr(m7) - ## Variance StdDev - ## block = pdLogChol(1) - ## (Intercept) 214.4716 14.64485 - ## gen = pdLogChol(1) - ## (Intercept) 109.6929 10.47344 - ## Residual 162.5591 12.74987 } + } } diff --git a/vignettes/agridat_data.Rmd b/vignettes/agridat_data.Rmd index 0bbd42fb..5209f135 100644 --- a/vignettes/agridat_data.Rmd +++ b/vignettes/agridat_data.Rmd @@ -5,7 +5,7 @@ output: rmarkdown::html_vignette: md_extensions: -autolink_bare_uris vignette: > - %\VignetteIndexEntry{Additional agricultural data} + %\VignetteIndexEntry{Additional sources of agricultural data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- diff --git a/vignettes/agridat_graphical_gems.Rmd b/vignettes/agridat_graphical_gems.Rmd index 18999411..171677cf 100644 --- a/vignettes/agridat_graphical_gems.Rmd +++ b/vignettes/agridat_graphical_gems.Rmd @@ -2,7 +2,7 @@ title: "Graphical Gems in the agridat Package" author: "Kevin Wright" output: - rmarkdown::html_vignette: + rmarkdown:html_vignette: md_extensions: -autolink_bare_uris bibliography: agridat.bib vignette: > diff --git a/vignettes/agridat_intro.Rmd b/vignettes/agridat_intro.Rmd index 19963dea..e2954278 100644 --- a/vignettes/agridat_intro.Rmd +++ b/vignettes/agridat_intro.Rmd @@ -1,5 +1,5 @@ --- -title: "Introduction to `agridat`" +title: "Introduction to agridat" author: "Kevin Wright" bibliography: agridat.bib output: