Single Site OFE
Data visualization
The data set is lasrosas.corn from the R-package agridat. The yield monitor data for a corn field (quintals/ha) was conducted by incorporating six nitrogen rate treatments in three replicated blocks comprising 18 strips across the field in Argentina in 2001. The raw data is converted to 93 rows \(\times\) 18 cols grid.
summary(gridded_df_ordered)## x y yield nitro
## Min. :-5137122 Min. :-6258764 Min. : 12.66 Min. : 0.0
## 1st Qu.:-5136910 1st Qu.:-6258727 1st Qu.: 50.51 1st Qu.: 39.0
## Median :-5136707 Median :-6258689 Median : 84.66 Median : 63.0
## Mean :-5136705 Mean :-6258689 Mean : 75.32 Mean : 64.9
## 3rd Qu.:-5136504 3rd Qu.:-6258650 3rd Qu.: 99.52 3rd Qu.: 99.8
## Max. :-5136293 Max. :-6258613 Max. :117.90 Max. :124.6
##
## topo bv rep nf StripNo RowId
## E :350 Min. : 91.74 R1:558 N0:279 Strip: 1: 93 Min. : 1
## HT:414 1st Qu.:168.10 R2:558 N1:279 Strip: 2: 93 1st Qu.:24
## LO:424 Median :173.63 R3:558 N2:279 Strip: 3: 93 Median :47
## W :486 Mean :174.67 N3:279 Strip: 4: 93 Mean :47
## 3rd Qu.:180.26 N4:279 Strip: 5: 93 3rd Qu.:70
## Max. :213.82 N5:279 Strip: 6: 93 Max. :93
## (Other) :1116
## row col cnitro nitro.sq cnitro.sq
## 1 : 18 1 : 93 Min. :-64.9 Min. : 0 Min. : 110.2
## 2 : 18 2 : 93 1st Qu.:-25.9 1st Qu.: 1521 1st Qu.: 204.5
## 3 : 18 3 : 93 Median : -1.9 Median : 4123 Median : 944.4
## 4 : 18 4 : 93 Mean : 0.0 Mean : 5875 Mean :1663.3
## 5 : 18 5 : 93 3rd Qu.: 34.9 3rd Qu.: 9960 3rd Qu.:3564.1
## 6 : 18 6 : 93 Max. : 59.7 Max. :15525 Max. :4212.0
## (Other):1566 (Other):1116
## gridobject gridId
## 1 : 6 1 : 1
## 2 : 6 2 : 1
## 3 : 6 3 : 1
## 4 : 6 4 : 1
## 5 : 6 5 : 1
## 6 : 6 6 : 1
## (Other):1638 (Other):1668ggplot(gridded_df_ordered) + geom_point(aes(row,col,colour=rep),size=2) + theme_bw()
ggplot(gridded_df_ordered) + geom_point(aes(row,col,colour=nf),size=2) + theme_bw()
ggplot(gridded_df_ordered) + geom_point(aes(row,col,colour=yield),size=2) + theme_bw()
ggplot(gridded_df_ordered) + geom_point(aes(row,col,colour=topo),size=2) + theme_bw()
ggplot(gridded_df_ordered, aes(x= yield,y=topo,fill=topo)) +
geom_density_ridges()## Picking joint bandwidth of 3.01
countplot <- gridded_df_ordered
countplot$x1<- as.numeric(countplot$row)
countplot$y1<- as.numeric(countplot$col)
coordinates(countplot) <- c("x1", "y1")
spplot(countplot, "yield")
countplot.ppp <- as(countplot["yield"], "ppp")
countplot.ppp.smooth <- Smooth(countplot.ppp)
plot(countplot.ppp.smooth, col=matlab.like2, ribsep=0.01,
ribside="bottom", main="")
contour(countplot.ppp.smooth,"marks.nitro", add=TRUE, lwd=1,
vfont=c("sans serif", "bold italic"), labcex=0.9)
Models
df <- read.csv("SinglesiteModel.csv")
kable(df)%>%kable_styling(bootstrap_options = c("striped", "hover"))| X | Fixed | Random | Residual | Log.likelihood | AIC | BIC |
|---|---|---|---|---|---|---|
| model 0 | 1 | rep+nf | — | -6266.143 | 12538.290 | 12554.550 |
| model 1 | 1 | rep+nf | id(row):id(col) | -6266.143 | 12538.290 | 12554.550 |
| model 2 | 1+lin(row)+lin(col) | rep+nf | ar1(row):ar1(col) | -3913.398 | 7836.796 | 7863.902 |
| model 3 | 1+lin(row) | rep+nf+spl(row) | ar1(row):ar1(col) | -3607.354 | 7226.708 | 7259.238 |
| model 4 | 1+lin(row) | rep+nf+spl(row)+outl | ar1(row):ar1(col) | -3441.974 | 6897.948 | 6935.900 |
| model 5 | 1+lin(row)+nf | rep+spl(row)+outl | ar1(row):ar1(col) | -3434.315 | 6880.630 | 6913.143 |
| model 6 | 1+lin(row)+nf+lin(row):nf | rep+spl(row)+outl | ar1(row):ar1(col) | -3450.218 | 6912.436 | 6944.931 |
| model 7 | 1+lin(row)+nf | rep+spl(row)+outl+nf:spl(row) | ar1(row):ar1(col) | -3430.177 | 6874.353 | 6912.284 |
| model 8 | 1+cnitro+cnitro.sq | rep | ar1(row):ar1(col) | -3922.920 | 7853.840 | 7875.525 |
| model 9 | 1+lin(row)+cnitro+cnitro.sq | rep+spl(row) | ar1(row):ar1(col) | -3611.849 | 7233.698 | 7260.801 |
| model 10 | 1+lin(row)+cnitro+cnitro.sq | rep+spl(row)+outl | ar1(row):ar1(col) | -3446.323 | 6904.646 | 6937.170 |
| model 11 | 1+lin(row)+cnitro+cnitro.sq | rep+spl(row)+outl+cnitro:spl(row)+cnitro.sq:spl(row) | ar1(row):ar1(col) | -3438.807 | 6893.614 | 6936.979 |
Fit treatment as random
#initial model
mod0 <- asreml(yield ~ 1,
random= ~ rep+nf,
data=gridded_df_ordered, maxiter=20)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:00 2020
## LogLik Sigma2 DF wall cpu
## 1 -6273.834 652.480 1673 09:37:00 0.0 (2 restrained)
## 2 -6267.399 654.012 1673 09:37:00 0.0 (1 restrained)
## 3 -6266.281 654.396 1673 09:37:00 0.0 (1 restrained)
## 4 -6266.152 654.576 1673 09:37:00 0.0 (1 restrained)
## 5 -6266.144 654.571 1673 09:37:00 0.0 (1 restrained)
## 6 -6266.143 654.571 1673 09:37:00 0.0asrmod(mod0)## $loglik
## [1] -6266.143
##
## $AIC
## [1] 12538.29
## attr(,"parameters")
## [1] 3
##
## $BIC
## [1] 12554.55
## attr(,"parameters")
## [1] 3plot(mod0)
mod1 <- asreml(yield ~ 1,
random= ~ rep+nf,
residual= ~ id(row):id(col),
data=gridded_df_ordered)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:02 2020
## LogLik Sigma2 DF wall cpu
## 1 -6273.834 652.480 1673 09:37:02 0.0 (2 restrained)
## 2 -6267.399 654.012 1673 09:37:02 0.0 (1 restrained)
## 3 -6266.281 654.396 1673 09:37:02 0.0 (1 restrained)
## 4 -6266.152 654.576 1673 09:37:02 0.0 (1 restrained)
## 5 -6266.144 654.571 1673 09:37:02 0.0 (1 restrained)
## 6 -6266.143 654.571 1673 09:37:02 0.0asrmod(mod1)## $loglik
## [1] -6266.143
##
## $AIC
## [1] 12538.29
## attr(,"parameters")
## [1] 3
##
## $BIC
## [1] 12554.55
## attr(,"parameters")
## [1] 3summary(mod1,param="gamma")$varcomp## gamma component std.error z.ratio bound %ch
## rep 1.011929e-07 6.623795e-05 2.293633e-06 28.8790582 B 0
## nf 4.658895e-03 3.049578e+00 3.413517e+00 0.8933832 P 0
## row:col!R 1.000000e+00 6.545712e+02 2.266595e+01 28.8790582 P 0#check the significance of the fixed effects
wald(mod1)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 4129124 6308.1 < 2.2e-16 ***
## residual (MS) 655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#run the diagnostics
plot(varioGram(mod1))
plot(mod1)
#fit the linear column and linear row effects model
mod2 <- asreml(yield ~ 1+lin(row)+lin(col),
random= ~ rep+nf,
residual= ~ar1(row):ar1(col),
data=gridded_df_ordered,maxiter=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:04 2020
## LogLik Sigma2 DF wall cpu
## 1 -5673.261 323.517 1671 09:37:04 0.0 (3 restrained)
## 2 -4615.786 126.868 1671 09:37:04 0.0 (2 restrained)
## 3 -3997.362 134.557 1671 09:37:04 0.0 (1 restrained)
## 4 -3963.913 161.589 1671 09:37:04 0.0 (1 restrained)
## 5 -3919.907 300.788 1671 09:37:04 0.0 (2 restrained)
## 6 -3913.614 345.522 1671 09:37:04 0.0 (1 restrained)
## 7 -3913.406 337.173 1671 09:37:04 0.0 (1 restrained)
## 8 -3913.398 341.281 1671 09:37:04 0.0
## 9 -3913.398 340.286 1671 09:37:04 0.0#diagnostics
asrmod(mod2)## $loglik
## [1] -3913.398
##
## $AIC
## [1] 7836.796
## attr(,"parameters")
## [1] 5
##
## $BIC
## [1] 7863.902
## attr(,"parameters")
## [1] 5# -3913.398
summary(mod2)$varcomp## component std.error z.ratio bound %ch
## rep 3.443454e-05 5.571003e-06 6.181031 B 0
## nf 8.883445e-05 1.437211e-05 6.181031 B 0
## row:col!R 3.402861e+02 5.505330e+01 6.181031 P 0
## row:col!row!cor 9.374610e-01 1.068600e-02 87.727923 U 0
## row:col!col!cor 2.639157e-01 2.502896e-02 10.544416 U 0wald(mod2)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 277682 816.03 < 2.2e-16 ***
## lin(row) 1 2981 8.76 0.003077 **
## lin(col) 1 37 0.11 0.743124
## residual (MS) 340
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod2))
plot(mod2)
# lin(col) is not significant, add spline row
mod3 <- asreml(yield~ 1+lin(row),
random= ~rep+nf+spl(row),
residual= ~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:06 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3736.861 30.5672 1672 09:37:06 0.1 (1 restrained)
## 2 -3671.982 29.3929 1672 09:37:06 0.1
## 3 -3624.633 29.7522 1672 09:37:07 0.1
## 4 -3609.884 31.0471 1672 09:37:07 0.1
## 5 -3607.429 32.0617 1672 09:37:07 0.1
## 6 -3607.356 32.2887 1672 09:37:07 0.1
## 7 -3607.354 32.3363 1672 09:37:07 0.1asrmod(mod3)#$loglik## $loglik
## [1] -3607.354
##
## $AIC
## [1] 7226.708
## attr(,"parameters")
## [1] 6
##
## $BIC
## [1] 7259.238
## attr(,"parameters")
## [1] 6# -3607.354
summary(mod3)$varcomp## component std.error z.ratio bound %ch
## rep 0.1837729 0.36505126 0.5034168 P 0.1
## nf 5.1187964 3.41221249 1.5001400 P 0.0
## spl(row) 9.8298245 2.96299498 3.3175299 P 0.0
## row:col!R 32.3362912 1.36943864 23.6128075 P 0.0
## row:col!row!cor 0.4415356 0.02081960 21.2076837 U 0.1
## row:col!col!cor 0.1486154 0.02568646 5.7857458 U 0.0wald(mod3)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 188245 5821.5 < 2.2e-16 ***
## lin(row) 1 91584 2832.2 < 2.2e-16 ***
## residual (MS) 32
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod3))
plot(mod3)
###############################################
# identify outliers
mod3.aom <- update(mod3,aom=T)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:09 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3607.354 32.3452 1672 09:37:09 0.1
## 2 -3607.354 32.3458 1672 09:37:09 0.1
## 3 -3607.354 32.3470 1672 09:37:09 0.1# gridded_df_ordered$aom <- residuals(mod3.aom)
gridded_df_ordered$aom <- mod3.aom$aom$R[,2]
ggplot(gridded_df_ordered) + geom_histogram(aes(aom),binwidth=1) +
ggtitle("Histogram of Standardised Residuals")+theme_bw()
# ggplot(gridded_df_ordered) + geom_boxplot(aes(topo,aom)) + theme_bw()
# standardised residual plots of model 3
mod3.1 <- mod3
mod3.1$residuals<- mod3.aom$aom$R[,2]
plot(mod3.1)
temp <- gridded_df_ordered
order(gridded_df_ordered$aom)[1:10]## [1] 74 1490 1116 1239 1369 1670 943 1473 1531 57wh <- which(abs(gridded_df_ordered$aom)>5)
resID <- c(74)
gridded_df_ordered[resID,]## x y yield nitro topo bv rep nf StripNo RowId row col
## 510 -5137046 -6258749 19.41 75.4 W 177.46 R1 N3 Strip: 2 5 5 2
## cnitro nitro.sq cnitro.sq gridobject gridId aom
## 510 10.5 5685.16 110.25 5 74 -17.228gridded_df_ordered$outl <- 0
gridded_df_ordered$outl[resID] <- 1
gridded_df_ordered$outl <- as.factor(gridded_df_ordered$outl)
# removing extreme outliers
gridded_df_ordered$yield[74] <- NA
mod3.2 <- asreml(yield~ 1+lin(row),
random= ~rep+nf+spl(row),
residual= ~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:11 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3628.468 26.9187 1671 09:37:11 0.1 (1 restrained)
## 2 -3528.873 25.1916 1671 09:37:11 0.1
## 3 -3460.926 25.5634 1671 09:37:11 0.1
## 4 -3440.607 27.1242 1671 09:37:11 0.1
## 5 -3437.212 28.4171 1671 09:37:11 0.1
## 6 -3437.106 28.7254 1671 09:37:11 0.1
## 7 -3437.102 28.7894 1671 09:37:11 0.1
## 8 -3437.102 28.8016 1671 09:37:11 0.1asrmod(mod3.2)#$loglik## $loglik
## [1] -3437.102
##
## $AIC
## [1] 6886.204
## attr(,"parameters")
## [1] 6
##
## $BIC
## [1] 6918.731
## attr(,"parameters")
## [1] 6summary(mod3.2)$varcomp## component std.error z.ratio bound %ch
## rep 0.1996678 0.39200006 0.5093566 P 0
## nf 5.2028157 3.47688797 1.4964002 P 0
## spl(row) 9.3337692 2.71290752 3.4405040 P 0
## row:col!R 28.8015509 1.29788494 22.1911435 P 0
## row:col!row!cor 0.5118910 0.02010364 25.4626049 U 0
## row:col!col!cor 0.1479415 0.02593579 5.7041445 U 0wald(mod3.2)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 163993 5693.9 < 2.2e-16 ***
## lin(row) 1 76586 2659.1 < 2.2e-16 ***
## residual (MS) 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod3.2))
plot(mod3.2)
mod3.2.aom <- update(mod3.2,aom=T)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:13 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3437.102 28.8039 1671 09:37:13 0.1
## 2 -3437.102 28.8041 1671 09:37:14 0.1
## 3 -3437.102 28.8044 1671 09:37:14 0.1gridded_df_ordered$aom2 <- mod3.2.aom$aom$R[,2]
ggplot(gridded_df_ordered) + geom_histogram(aes(aom2),binwidth=1) +
ggtitle("Histogram of Standardised Residuals")+theme_bw()
order(gridded_df_ordered$aom2)[1:10]## [1] 1490 1116 1239 1369 1670 943 1473 719 1131 1531resID <- c(74,which(abs(gridded_df_ordered$aom2)>2))
gridded_df_ordered[resID,]## x y yield nitro topo bv rep nf StripNo RowId row col
## 510 -5137046 -6258749 NA 75.4 W 177.46 R1 N3 Strip: 2 5 5 2
## 110 -5137080 -6258749 87.73 75.4 W 166.49 R1 N3 Strip: 2 1 1 2
## 1101 -5137122 -6258613 86.66 39.0 W 142.55 R3 N1 Strip: 18 1 1 18
## 270 -5137105 -6258647 100.52 75.4 W 139.78 R3 N3 Strip: 14 2 2 14
## 2101 -5137113 -6258613 102.16 39.0 W 131.17 R3 N1 Strip: 18 2 2 18
## 412 -5137063 -6258741 91.04 99.8 W 174.33 R1 N4 Strip: 3 4 4 3
## 460 -5137088 -6258656 110.90 124.6 W 96.39 R3 N5 Strip: 13 4 4 13
## 570 -5137080 -6258647 80.05 75.4 W 171.65 R3 N3 Strip: 14 5 5 14
## 580 -5137079 -6258639 85.22 99.8 W 172.75 R3 N4 Strip: 15 5 5 15
## 5101 -5137088 -6258614 99.69 39.0 W 171.22 R3 N1 Strip: 18 5 5 18
## 610 -5137038 -6258749 108.03 75.4 W 173.70 R1 N3 Strip: 2 6 6 2
## 650 -5137071 -6258639 108.45 99.8 W 174.05 R3 N4 Strip: 15 6 6 15
## 812 -5137029 -6258741 100.54 99.8 W 177.22 R1 N4 Strip: 3 8 8 3
## 9 -5137013 -6258758 100.44 124.6 W 179.49 R1 N5 Strip: 1 9 9 1
## 920 -5137038 -6258690 76.03 99.8 W 176.80 R2 N4 Strip: 9 9 9 9
## 1112 -5137021 -6258690 78.20 99.8 W 175.24 R2 N4 Strip: 9 11 11 9
## 121 -5136987 -6258750 81.90 75.4 W 183.94 R1 N3 Strip: 2 12 12 2
## 128 -5137004 -6258699 49.13 75.4 W 185.62 R2 N3 Strip: 8 12 12 8
## 13 -5136979 -6258758 82.70 124.6 W 180.27 R1 N5 Strip: 1 13 13 1
## 23 -5136894 -6258759 62.98 124.6 W 181.74 R1 N5 Strip: 1 23 23 1
## 259 -5136902 -6258683 36.08 0.0 W 177.30 R2 N0 Strip: 10 25 25 10
## 2615 -5136910 -6258632 36.09 0.0 W 181.43 R3 N0 Strip: 16 26 26 16
## 3315 -5136851 -6258632 31.38 0.0 HT 171.74 R3 N0 Strip: 16 33 33 16
## 3614 -5136817 -6258641 37.03 99.8 HT 188.19 R3 N4 Strip: 15 35 35 15
## 3917 -5136800 -6258616 52.90 39.0 HT 202.74 R3 N1 Strip: 18 38 38 18
## 4013 -5136783 -6258650 34.62 75.4 HT 198.83 R3 N3 Strip: 14 39 39 14
## 4117 -5136775 -6258650 43.70 75.4 HT 203.71 R3 N3 Strip: 14 40 40 14
## 4118 -5136775 -6258641 32.05 99.8 HT 204.97 R3 N4 Strip: 15 40 40 15
## 4120 -5136783 -6258624 31.79 50.6 HT 203.68 R3 N2 Strip: 17 40 40 17
## 4218 -5136775 -6258625 47.02 50.6 HT 202.14 R3 N2 Strip: 17 41 41 17
## 423 -5136742 -6258735 12.66 0.0 HT 193.16 R1 N0 Strip: 4 42 42 4
## 442 -5136725 -6258743 32.90 99.8 HT 187.79 R1 N4 Strip: 3 44 44 3
## 4516 -5136749 -6258625 50.57 50.6 HT 182.88 R3 N2 Strip: 17 44 44 17
## 4713 -5136732 -6258616 62.88 39.0 HT 179.37 R3 N1 Strip: 18 46 46 18
## 5218 -5136690 -6258625 68.11 50.6 HT 182.17 R3 N2 Strip: 17 50 50 17
## 5212 -5136673 -6258676 47.52 50.6 E 177.43 R2 N2 Strip: 11 51 51 11
## 531 -5136640 -6258753 75.83 75.4 E 177.89 R1 N3 Strip: 2 52 52 2
## 546 -5136648 -6258710 58.28 124.6 E 176.14 R2 N5 Strip: 7 53 53 7
## 547 -5136648 -6258702 74.35 75.4 E 175.50 R2 N3 Strip: 8 53 53 8
## 548 -5136657 -6258693 63.75 99.8 E 178.95 R2 N4 Strip: 9 53 53 9
## 556 -5136640 -6258710 87.78 124.6 E 177.27 R2 N5 Strip: 7 54 54 7
## 5415 -5136648 -6258634 69.74 0.0 E 173.77 R3 N0 Strip: 16 54 54 16
## 5614 -5136648 -6258642 65.65 99.8 E 177.05 R3 N4 Strip: 15 55 55 15
## 5717 -5136648 -6258617 81.37 39.0 E 173.73 R3 N1 Strip: 18 55 55 18
## 5817 -5136639 -6258617 101.34 39.0 E 174.48 R3 N1 Strip: 18 56 56 18
## 583 -5136606 -6258736 103.97 0.0 E 170.44 R1 N0 Strip: 4 57 57 4
## 5915 -5136605 -6258634 91.11 0.0 E 171.43 R3 N0 Strip: 16 59 59 16
## 6118 -5136605 -6258643 98.81 99.8 E 171.22 R3 N4 Strip: 15 60 60 15
## 6015 -5136597 -6258634 102.16 0.0 E 171.42 R3 N0 Strip: 16 60 60 16
## 635 -5136572 -6258719 105.73 39.0 E 176.47 R1 N1 Strip: 6 61 61 6
## 645 -5136564 -6258719 97.94 39.0 E 177.19 R1 N1 Strip: 6 62 62 6
## 6414 -5136580 -6258643 111.74 99.8 E 175.66 R3 N4 Strip: 15 62 62 15
## 6416 -5136588 -6258626 109.24 50.6 E 164.57 R3 N2 Strip: 17 62 62 17
## 6417 -5136588 -6258618 81.00 39.0 E 163.03 R3 N1 Strip: 18 62 62 18
## 655 -5136555 -6258719 106.84 39.0 E 169.45 R1 N1 Strip: 6 63 63 6
## 6514 -5136572 -6258643 95.77 99.8 E 167.74 R3 N4 Strip: 15 63 63 15
## 6516 -5136580 -6258626 92.21 50.6 E 168.51 R3 N2 Strip: 17 63 63 17
## 6811 -5136538 -6258669 110.54 39.0 E 161.69 R2 N1 Strip: 12 66 66 12
## 6916 -5136546 -6258626 96.32 50.6 E 167.79 R3 N2 Strip: 17 67 67 17
## 708 -5136521 -6258694 98.09 99.8 E 169.45 R2 N4 Strip: 9 68 68 9
## 7016 -5136538 -6258626 112.85 50.6 E 170.98 R3 N2 Strip: 17 68 68 17
## 7112 -5136513 -6258694 111.87 99.8 E 171.61 R2 N4 Strip: 9 69 69 9
## 7118 -5136521 -6258643 90.45 99.8 E 169.94 R3 N4 Strip: 15 69 69 15
## 7313 -5136504 -6258652 94.33 75.4 E 165.49 R3 N3 Strip: 14 71 71 14
## 734 -5136479 -6258728 93.06 50.6 LO 169.12 R1 N2 Strip: 5 72 72 5
## 7414 -5136495 -6258644 98.88 99.8 E 161.76 R3 N4 Strip: 15 72 72 15
## 744 -5136471 -6258729 108.01 50.6 LO 162.02 R1 N2 Strip: 5 73 73 5
## 7614 -5136478 -6258644 117.90 99.8 LO 162.17 R3 N4 Strip: 15 74 74 15
## 76 -5136446 -6258763 95.93 124.6 LO 169.89 R1 N5 Strip: 1 75 75 1
## 77 -5136437 -6258763 116.64 124.6 LO 159.75 R1 N5 Strip: 1 76 76 1
## 78 -5136429 -6258763 97.56 124.6 LO 163.36 R1 N5 Strip: 1 77 77 1
## 7914 -5136453 -6258644 112.94 99.8 LO 172.52 R3 N4 Strip: 15 77 77 15
## 8113 -5136428 -6258687 107.59 0.0 LO 169.97 R2 N0 Strip: 10 79 79 10
## 8118 -5136436 -6258644 87.43 99.8 LO 163.91 R3 N4 Strip: 15 79 79 15
## 813 -5136412 -6258746 112.83 99.8 LO 170.25 R1 N4 Strip: 3 80 80 3
## 8211 -5136420 -6258687 94.65 0.0 LO 165.06 R2 N0 Strip: 10 80 80 10
## 8417 -5136419 -6258619 107.22 39.0 LO 164.82 R3 N1 Strip: 18 81 81 18
## 83 -5136386 -6258763 110.56 124.6 LO 163.70 R1 N5 Strip: 1 82 82 1
## 8413 -5136411 -6258653 100.51 75.4 LO 153.08 R3 N3 Strip: 14 82 82 14
## 8414 -5136411 -6258644 92.60 99.8 LO 161.18 R3 N4 Strip: 15 82 82 15
## 8513 -5136402 -6258653 69.21 75.4 LO 163.36 R3 N3 Strip: 14 83 83 14
## 8514 -5136402 -6258644 106.21 99.8 LO 167.98 R3 N4 Strip: 15 83 83 15
## 8315 -5136402 -6258636 102.93 0.0 LO 165.34 R3 N0 Strip: 16 83 83 16
## 8612 -5136394 -6258661 102.66 124.6 LO 165.70 R3 N5 Strip: 13 84 84 13
## 8613 -5136394 -6258653 105.25 75.4 LO 164.87 R3 N3 Strip: 14 84 84 14
## 8415 -5136394 -6258636 86.67 0.0 LO 161.98 R3 N0 Strip: 16 84 84 16
## 862 -5136369 -6258746 109.52 99.8 LO 161.32 R1 N4 Strip: 3 85 85 3
## 8712 -5136385 -6258662 117.19 124.6 LO 165.81 R3 N5 Strip: 13 85 85 13
## 87 -5136352 -6258763 85.64 124.6 LO 169.63 R1 N5 Strip: 1 86 86 1
## 871 -5136352 -6258755 102.12 75.4 LO 177.52 R1 N3 Strip: 2 86 86 2
## 872 -5136361 -6258746 94.48 99.8 LO 171.21 R1 N4 Strip: 3 86 86 3
## 9214 -5136343 -6258653 109.57 75.4 LO 164.81 R3 N3 Strip: 14 90 90 14
## 9015 -5136343 -6258636 85.13 0.0 LO 167.75 R3 N0 Strip: 16 90 90 16
## 9316 -5136343 -6258628 100.36 50.6 LO 166.01 R3 N2 Strip: 17 90 90 17
## 9413 -5136326 -6258654 104.21 75.4 LO 173.38 R3 N3 Strip: 14 92 92 14
## 9414 -5136326 -6258645 107.01 99.8 LO 170.19 R3 N4 Strip: 15 92 92 15
## 959 -5136318 -6258654 85.05 75.4 LO 167.84 R3 N3 Strip: 14 93 93 14
## cnitro nitro.sq cnitro.sq gridobject gridId aom outl aom2
## 510 10.5 5685.16 110.25 5 74 -17.228003 1 NA
## 110 10.5 5685.16 110.25 1 2 -1.810830 0 -2.042881
## 1101 -25.9 1521.00 670.81 187 18 -2.074970 0 -2.364591
## 270 10.5 5685.16 110.25 188 32 1.906677 0 2.251556
## 2101 -25.9 1521.00 670.81 188 36 2.386740 0 2.832189
## 412 34.9 9960.04 1218.01 4 57 -2.598582 0 -2.003388
## 460 59.7 15525.16 3564.09 190 67 2.964006 0 3.252966
## 570 10.5 5685.16 110.25 191 86 -2.071350 0 -2.365493
## 580 34.9 9960.04 1218.01 191 87 -1.774945 0 -2.112016
## 5101 -25.9 1521.00 670.81 191 90 2.328933 0 2.637118
## 610 10.5 5685.16 110.25 6 92 8.586755 0 2.611539
## 650 34.9 9960.04 1218.01 192 105 3.148979 0 3.561319
## 812 34.9 9960.04 1218.01 8 129 1.981116 0 2.234895
## 9 59.7 15525.16 3564.09 9 145 1.914733 0 2.045011
## 920 34.9 9960.04 1218.01 102 153 -1.770748 0 -2.023645
## 1112 34.9 9960.04 1218.01 104 189 1.763364 0 2.018685
## 121 10.5 5685.16 110.25 12 200 2.374815 0 2.631044
## 128 10.5 5685.16 110.25 105 206 -2.209220 0 -2.415099
## 13 59.7 15525.16 3564.09 13 217 2.432004 0 2.556345
## 23 59.7 15525.16 3564.09 23 397 1.787428 0 2.033787
## 259 -64.9 0.00 4212.01 118 442 -1.852269 0 -2.063370
## 2615 -64.9 0.00 4212.01 212 466 -2.089166 0 -2.315745
## 3315 -64.9 0.00 4212.01 219 592 -1.935853 0 -2.125031
## 3614 34.9 9960.04 1218.01 221 627 -2.301760 0 -2.603531
## 3917 -25.9 1521.00 670.81 224 684 2.200645 0 2.429073
## 4013 10.5 5685.16 110.25 225 698 -1.889680 0 -2.295096
## 4117 10.5 5685.16 110.25 226 716 1.823544 0 2.056192
## 4118 34.9 9960.04 1218.01 226 717 -1.865773 0 -2.150357
## 4120 -14.3 2560.36 204.49 226 719 -2.561293 0 -3.077060
## 4218 -14.3 2560.36 204.49 227 737 3.090783 0 3.575272
## 423 -64.9 0.00 4212.01 42 742 -2.344173 0 -2.442339
## 442 34.9 9960.04 1218.01 44 777 -1.795709 0 -2.028226
## 4516 -14.3 2560.36 204.49 230 791 1.921436 0 2.151343
## 4713 -25.9 1521.00 670.81 232 828 2.241009 0 2.415864
## 5218 -14.3 2560.36 204.49 236 899 1.981692 0 2.176197
## 5212 -14.3 2560.36 204.49 144 911 -2.063355 0 -2.373900
## 531 10.5 5685.16 110.25 52 920 1.933792 0 2.130680
## 546 59.7 15525.16 3564.09 146 943 -2.868999 0 -3.267559
## 547 10.5 5685.16 110.25 146 944 1.753662 0 2.001848
## 548 34.9 9960.04 1218.01 146 945 -2.124362 0 -2.448849
## 556 59.7 15525.16 3564.09 147 961 2.705213 0 3.176615
## 5415 -64.9 0.00 4212.01 240 970 -1.801181 0 -2.188551
## 5614 34.9 9960.04 1218.01 241 987 -2.175815 0 -2.228924
## 5717 -25.9 1521.00 670.81 241 990 -2.083140 0 -2.502616
## 5817 -25.9 1521.00 670.81 242 1008 2.268576 0 2.575574
## 583 -64.9 0.00 4212.01 57 1012 1.896104 0 2.060612
## 5915 -64.9 0.00 4212.01 245 1060 -1.882338 0 -2.188605
## 6118 34.9 9960.04 1218.01 246 1077 -1.733118 0 -2.016476
## 6015 -64.9 0.00 4212.01 246 1078 2.232872 0 2.659047
## 635 -25.9 1521.00 670.81 61 1086 1.692664 0 2.002704
## 645 -25.9 1521.00 670.81 62 1104 -1.723047 0 -2.054753
## 6414 34.9 9960.04 1218.01 248 1113 1.863626 0 2.148763
## 6416 -14.3 2560.36 204.49 248 1115 1.886120 0 2.103782
## 6417 -25.9 1521.00 670.81 248 1116 -3.462416 0 -3.804494
## 655 -25.9 1521.00 670.81 63 1122 1.764086 0 2.083405
## 6514 34.9 9960.04 1218.01 249 1131 -2.527782 0 -2.967026
## 6516 -14.3 2560.36 204.49 249 1133 -2.164896 0 -2.500754
## 6811 -25.9 1521.00 670.81 159 1182 1.909261 0 2.140685
## 6916 -14.3 2560.36 204.49 253 1205 -2.141243 0 -2.525669
## 708 34.9 9960.04 1218.01 161 1215 -1.805094 0 -2.112972
## 7016 -14.3 2560.36 204.49 254 1223 2.765228 0 3.202909
## 7112 34.9 9960.04 1218.01 162 1233 2.021065 0 2.331553
## 7118 34.9 9960.04 1218.01 255 1239 -3.225017 0 -3.680258
## 7313 10.5 5685.16 110.25 257 1274 -2.387866 0 -2.759634
## 734 -14.3 2560.36 204.49 72 1283 -2.505848 0 -2.917365
## 7414 34.9 9960.04 1218.01 258 1293 -2.050650 0 -2.386792
## 744 -14.3 2560.36 204.49 73 1301 1.928096 0 2.254462
## 7614 34.9 9960.04 1218.01 260 1329 2.334467 0 2.602642
## 76 59.7 15525.16 3564.09 75 1333 -2.316430 0 -2.687370
## 77 59.7 15525.16 3564.09 76 1351 3.710383 0 4.346406
## 78 59.7 15525.16 3564.09 77 1369 -3.022961 0 -3.621419
## 7914 34.9 9960.04 1218.01 263 1383 2.880488 0 3.350964
## 8113 -64.9 0.00 4212.01 172 1414 2.063600 0 2.345742
## 8118 34.9 9960.04 1218.01 265 1419 -2.180484 0 -2.369311
## 813 34.9 9960.04 1218.01 80 1425 3.092094 0 3.596567
## 8211 -64.9 0.00 4212.01 173 1432 -2.005950 0 -2.420092
## 8417 -25.9 1521.00 670.81 267 1458 1.795682 0 2.103211
## 83 59.7 15525.16 3564.09 82 1459 2.009552 0 2.317294
## 8413 10.5 5685.16 110.25 268 1472 3.234469 0 4.033093
## 8414 34.9 9960.04 1218.01 268 1473 -2.700483 0 -3.117826
## 8513 10.5 5685.16 110.25 269 1490 -7.622807 0 -8.716190
## 8514 34.9 9960.04 1218.01 269 1491 2.846929 0 3.368931
## 8315 -64.9 0.00 4212.01 269 1492 1.824425 0 2.138103
## 8612 59.7 15525.16 3564.09 270 1507 -1.784264 0 -2.209579
## 8613 10.5 5685.16 110.25 270 1508 3.757767 0 4.556601
## 8415 -64.9 0.00 4212.01 270 1510 -2.387898 0 -2.780062
## 862 34.9 9960.04 1218.01 85 1515 2.280352 0 2.625644
## 8712 59.7 15525.16 3564.09 271 1525 2.323150 0 2.534360
## 87 59.7 15525.16 3564.09 86 1531 -2.662678 0 -2.950584
## 871 10.5 5685.16 110.25 86 1532 1.750922 0 2.023950
## 872 34.9 9960.04 1218.01 86 1533 -1.775027 0 -2.059063
## 9214 10.5 5685.16 110.25 276 1616 2.418454 0 2.729257
## 9015 -64.9 0.00 4212.01 276 1618 -2.223945 0 -2.570816
## 9316 -14.3 2560.36 204.49 276 1619 1.943662 0 2.287110
## 9413 10.5 5685.16 110.25 278 1652 2.184576 0 2.594662
## 9414 34.9 9960.04 1218.01 278 1653 1.883662 0 2.186259
## 959 10.5 5685.16 110.25 279 1670 -3.021050 0 -3.436206gridded_df_ordered$outl <- 0
gridded_df_ordered$outl[resID] <- 1
gridded_df_ordered$outl <- as.factor(gridded_df_ordered$outl)
### re-fit without outliers
mod4 <- asreml(yield~ 1+lin(row),
random= ~rep+nf+spl(row)+outl,
residual= ~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:15 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3629.874 26.9174 1671 09:37:15 0.1 (2 restrained)
## 2 -3529.028 25.1879 1671 09:37:15 0.1 (1 restrained)
## 3 -3460.859 25.5663 1671 09:37:15 0.1 (1 restrained)
## 4 -3440.593 27.1274 1671 09:37:15 0.1 (1 restrained)
## 5 -3437.212 28.4179 1671 09:37:15 0.1 (1 restrained)
## 6 -3437.106 28.7256 1671 09:37:15 0.1
## 7 -3437.102 28.7894 1671 09:37:15 0.1
## 8 -3437.102 28.8016 1671 09:37:15 0.1asrmod(mod4)# $loglik## $loglik
## [1] -3437.102
##
## $AIC
## [1] 6888.204
## attr(,"parameters")
## [1] 7
##
## $BIC
## [1] 6926.152
## attr(,"parameters")
## [1] 7summary(mod4)$varcomp## component std.error z.ratio bound %ch
## outl 2.914513e-06 1.313368e-07 22.1911349 B 0
## rep 1.996674e-01 3.919996e-01 0.5093561 P 0
## nf 5.202816e+00 3.476888e+00 1.4964004 P 0
## spl(row) 9.333772e+00 2.712907e+00 3.4405054 P 0
## row:col!R 2.880156e+01 1.297886e+00 22.1911349 P 0
## row:col!row!cor 5.118911e-01 2.010364e-02 25.4626056 U 0
## row:col!col!cor 1.479415e-01 2.593579e-02 5.7041444 U 0wald(mod4)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 163993 5693.9 < 2.2e-16 ***
## lin(row) 1 76586 2659.1 < 2.2e-16 ***
## residual (MS) 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod4))
plot(mod4)
mod4.aom <- update(mod4,aom=T)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:17 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3437.102 28.8039 1671 09:37:17 0.1
## 2 -3437.102 28.8041 1671 09:37:17 0.1
## 3 -3437.102 28.8044 1671 09:37:17 0.1# standardised residual plots of model 4
mod4.1 <- mod4
mod4.1$residuals<- mod4.aom$aom$R[,2]
plot(mod4.1)
## BLUP
pv.blup4 <- predict(mod4,classify="nf")## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:19 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3437.102 28.8039 1671 09:37:19 0.1
## 2 -3437.102 28.8041 1671 09:37:19 0.1
## 3 -3437.102 28.8044 1671 09:37:19 0.1pv.blup4$pvals##
## Notes:
## [0;34m- The predictions are obtained by averaging across the hypertable
## calculated from model terms constructed solely from factors in
## the averaging and classify sets.[0m
## [0;34m- Use 'average' to move ignored factors into the averaging set.[0m
## [0;34m- lin(row) evaluated at average value of 47.000000[0m
## [0;34m- The ignored set: rep,outl[0m
## [0;34m
## [0m
## nf predicted.value std.error status
## 1 N0 43.75732 1.230931 Estimable
## 2 N1 47.02025 1.230343 Estimable
## 3 N2 47.67621 1.230902 Estimable
## 4 N3 48.65409 1.230951 Estimable
## 5 N4 49.45452 1.230932 Estimable
## 6 N5 49.85577 1.230343 Estimablelrt(mod3,mod4)## Likelihood ratio test(s) assuming nested random models.
## (See Self & Liang, 1987)
##
## df LR-statistic Pr(Chisq)
## mod4/mod3 1 340.5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fit treatment as fixed
mod5 <- asreml(yield~ 1+lin(row)+nf,
random= ~rep+spl(row)+outl,
residual=~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:24 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3621.512 26.8404 1666 09:37:24 0.1 (2 restrained)
## 2 -3521.236 25.1403 1666 09:37:24 0.1 (1 restrained)
## 3 -3453.252 25.5384 1666 09:37:24 0.1 (1 restrained)
## 4 -3432.959 27.1156 1666 09:37:24 0.1 (1 restrained)
## 5 -3429.555 28.4169 1666 09:37:24 0.1 (1 restrained)
## 6 -3429.448 28.7258 1666 09:37:24 0.1
## 7 -3429.444 28.7899 1666 09:37:24 0.1
## 8 -3429.444 28.8022 1666 09:37:24 0.1mod5 <- update(mod5)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:25 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3429.444 28.8046 1666 09:37:25 0.1
## 2 -3429.444 28.8047 1666 09:37:26 0.1asrmod(mod5)## $loglik
## [1] -3429.444
##
## $AIC
## [1] 6870.888
## attr(,"parameters")
## [1] 6
##
## $BIC
## [1] 6903.397
## attr(,"parameters")
## [1] 6summary(mod5)$varcomp## component std.error z.ratio bound %ch
## outl 2.914831e-06 1.313756e-07 22.1870037 B 0
## rep 1.986928e-01 3.910214e-01 0.5081378 P 0
## spl(row) 9.334449e+00 2.712630e+00 3.4411066 P 0
## row:col!R 2.880470e+01 1.298269e+00 22.1870037 P 0
## row:col!row!cor 5.119366e-01 2.010437e-02 25.4639431 U 0
## row:col!col!cor 1.478615e-01 2.593447e-02 5.7013516 U 0wald(mod5)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1216030 42216 < 2.2e-16 ***
## lin(row) 1 76582 2659 < 2.2e-16 ***
## nf 5 2759 96 < 2.2e-16 ***
## residual (MS) 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod5))
plot(mod5)
mod6 <- asreml(yield~ 1+lin(row)+nf+lin(row):nf,
random= ~rep+spl(row)+outl,
residual=~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:28 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3635.245 26.6853 1661 09:37:28 0.1 (2 restrained)
## 2 -3536.439 25.0406 1661 09:37:28 0.1 (1 restrained)
## 3 -3469.082 25.4564 1661 09:37:28 0.1 (1 restrained)
## 4 -3448.874 27.0394 1661 09:37:28 0.1 (1 restrained)
## 5 -3445.461 28.3482 1661 09:37:28 0.1 (1 restrained)
## 6 -3445.351 28.6624 1661 09:37:28 0.1
## 7 -3445.347 28.7289 1661 09:37:28 0.1
## 8 -3445.347 28.7418 1661 09:37:28 0.1mod6 <- update(mod6)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:29 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3445.347 28.7444 1661 09:37:29 0.1
## 2 -3445.347 28.7445 1661 09:37:29 0.1asrmod(mod6)## $loglik
## [1] -3445.347
##
## $AIC
## [1] 6902.694
## attr(,"parameters")
## [1] 6
##
## $BIC
## [1] 6935.185
## attr(,"parameters")
## [1] 6# -3450.218
summary(mod6)$varcomp## component std.error z.ratio bound %ch
## outl 2.908741e-06 1.317731e-07 22.0738570 B 0
## rep 1.987786e-01 3.905289e-01 0.5089984 P 0
## spl(row) 9.329717e+00 2.713675e+00 3.4380380 P 0
## row:col!R 2.874452e+01 1.302197e+00 22.0738570 P 0
## row:col!row!cor 5.100336e-01 2.035738e-02 25.0539977 U 0
## row:col!col!cor 1.511061e-01 2.604523e-02 5.8016782 U 0wald(mod6)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1213847 42229 <2e-16 ***
## lin(row) 1 76464 2660 <2e-16 ***
## nf 5 2781 97 <2e-16 ***
## lin(row):nf 5 167 6 0.3252
## residual (MS) 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod6))
plot(mod6)
mod7 <- asreml(yield~ 1+lin(row)+nf,
random= ~rep+outl+spl(row)+nf:spl(row),
residual=~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:31 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3650.042 24.2461 1666 09:37:32 0.3 (3 restrained)
## 2 -3526.211 23.8951 1666 09:37:32 0.3 (2 restrained)
## 3 -3449.071 24.6911 1666 09:37:32 0.3 (2 restrained)
## 4 -3430.819 26.4861 1666 09:37:32 0.3 (1 restrained)
## 5 -3425.773 27.6115 1666 09:37:33 0.3 (1 restrained)
## 6 -3425.339 27.8155 1666 09:37:33 0.3
## 7 -3425.309 27.8430 1666 09:37:33 0.3
## 8 -3425.309 27.8493 1666 09:37:33 0.3mod7 <- update(mod7)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:35 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3425.309 27.8508 1666 09:37:35 0.3
## 2 -3425.309 27.8509 1666 09:37:35 0.3asrmod(mod7)## $loglik
## [1] -3425.309
##
## $AIC
## [1] 6864.618
## attr(,"parameters")
## [1] 7
##
## $BIC
## [1] 6902.545
## attr(,"parameters")
## [1] 7summary(mod7)$varcomp## component std.error z.ratio bound %ch
## outl 2.818312e-06 1.268550e-07 22.2168025 B 0
## rep 2.117725e-01 3.908947e-01 0.5417634 P 0
## spl(row) 9.389952e+00 2.723416e+00 3.4478580 P 0
## nf:spl(row) 3.043402e-03 2.337963e-03 1.3017327 P 0
## row:col!R 2.785089e+01 1.253596e+00 22.2168025 P 0
## row:col!row!cor 4.933503e-01 2.094974e-02 23.5492360 U 0
## row:col!col!cor 1.565765e-01 2.605421e-02 6.0096423 U 0wald(mod7)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1176085 42228 < 2.2e-16 ***
## lin(row) 1 78843 2831 < 2.2e-16 ***
## nf 5 2926 105 < 2.2e-16 ***
## residual (MS) 28
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod7))
plot(mod7)
pv.blue7 <- predict(mod7,classify="nf")## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:37 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3425.309 27.8511 1666 09:37:37 0.3
## 2 -3425.309 27.8511 1666 09:37:38 0.3
## 3 -3425.309 27.8511 1666 09:37:38 0.3pv.blue7$pvals##
## Notes:
## [0;34m- The predictions are obtained by averaging across the hypertable
## calculated from model terms constructed solely from factors in
## the averaging and classify sets.[0m
## [0;34m- Use 'average' to move ignored factors into the averaging set.[0m
## [0;34m- lin(row) evaluated at average value of 47.000000[0m
## [0;34m- The ignored set: rep,outl[0m
## [0;34m
## [0m
## nf predicted.value std.error status
## 1 N0 42.28202 1.327859 Estimable
## 2 N1 47.89010 1.325923 Estimable
## 3 N2 48.60999 1.327754 Estimable
## 4 N3 49.20231 1.327757 Estimable
## 5 N4 49.47227 1.327859 Estimable
## 6 N5 49.02294 1.325924 Estimablelrt(mod5,mod7)## Likelihood ratio test(s) assuming nested random models.
## (See Self & Liang, 1987)
##
## df LR-statistic Pr(Chisq)
## mod7/mod5 1 8.2694 0.002016 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Compare the BLUEs, BLUPs and raw means
BLUE is from model 7 and BLUP is from model 4.
library(plyr)## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following object is masked from 'package:maps':
##
## ozone## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarizepv.blup.df=pv.blup4$pvals
pv.blue.df=pv.blue7$pvals
raw.mean=ddply(gridded_df_ordered,~nf,
summarise,mean=mean(yield, na.rm = TRUE),
sem=round(sd(yield, na.rm = TRUE)/sqrt(length(yield[!is.na(yield)])), 2))
pv.blup.df$status=NULL
pv.blue.df$status=NULL
pv.blup.df$type="BLUP"
pv.blue.df$type="BLUE"
raw.mean$type="RAW"
names(pv.blup.df)=c("Variety", "mean", "std.error", "type")
names(pv.blue.df)=c("Variety", "mean", "std.error", "type")
names(raw.mean)=c("Variety", "mean", "std.error", "type")
lup.pv1=cbind(raw.mean, pv.blup.df, pv.blue.df)
names(lup.pv1)=c("Variety","mean","std.error1", "type1","Variety2","BLUP","std.error2", "type2","Variety3",
"BLUE","std.error3", "type3" )
lup.pv2=rbind(raw.mean, pv.blup.df, pv.blue.df)
lup.pv2=lup.pv2[order(lup.pv2$type,lup.pv2$Variety, lup.pv2$mean),]
xyplot(mean~as.numeric(Variety),type=c('p'),data=lup.pv2,lwd=2,
main="Comparison", group=type,
xlab="Nitrogen Level",
ylab="Means",auto.key=list(columns=3),as.table=TRUE,
scales=list(x=list(at=seq(1,54), rot=90, cex=0.7, labels=levels(lup.pv2$Variety))),
par.settings=simpleTheme(col=c("red", "gold", "purple"),
pch=c(8, 6, 5),lty=(1:3), cex=0.5))
xyplot(mean~as.numeric(Variety),type=c('p', 'h'),data=lup.pv2,lwd=2,main="Comparison", group=type,
xlab="Nitrogen Level",
ylab="Means",auto.key=list(columns=3),as.table=TRUE,
scales=list(x=list(at=seq(1,54), rot=90, cex=0.7, labels=levels(lup.pv2$Variety))),
par.settings=simpleTheme(col=c("red", "gold", "purple"),
pch=c(8, 6, 5),lty=(1:3), cex=0.5))
lup.pv1$Variety=factor(lup.pv1$Variety)
library(lattice)
xyplot(BLUP+BLUE~mean,type=c('p', 'r'),data=lup.pv1,lwd=2,main="Comparison",
xlab="Observed means",
ylab="Predicted means",auto.key=list(columns=2),as.table=TRUE,
par.settings=simpleTheme(col=c("red", "blue"),
pch=c(8, 6),lty=(1:2), cex=0.5))
library(ggplot2)
pd <- position_dodge(width = 0.5)
ggplot(data=lup.pv2, aes(x = Variety, y = mean, group = type, colour=type)) +
geom_line(aes(linetype = type), position = pd, lwd=0.8) + scale_color_manual(values=c("red", "blue", "gold")) +
geom_point(size = 2, position = pd, aes(shape=type)) + theme(legend.position="top",axis.text.x=element_text(angle=90, hjust=1)) + theme_bw()
ggplot(data =lup.pv2, aes(x = type, y = mean, fill=type)) +
geom_boxplot(aes(fill = type), width = 0.85, position=position_dodge(width=0.85)) +
theme(legend.position="right") +
labs(y="Yield",x="") +
# scale_fill_manual(values = hmcol) +
stat_summary(fun.y="mean", colour="darkblue", geom="point",
position=position_dodge(width=0.85), aes(fill=type),
shape=18, size=2.5)+theme_bw()
Fit quadratic form \(\beta_0+\beta_1N+\beta_2N^2\)
mod8 <- asreml(yield~ 1+cnitro+cnitro.sq,
random= ~rep,
residual=~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:42 2020
## LogLik Sigma2 DF wall cpu
## 1 -5941.452 447.779 1670 09:37:42 0.0 (3 restrained)
## 2 -4781.814 154.203 1670 09:37:42 0.0 (2 restrained)
## 3 -3858.124 134.311 1670 09:37:42 0.0 (1 restrained)
## 4 -3816.302 161.165 1670 09:37:42 0.0 (2 restrained)
## 5 -3755.051 314.682 1670 09:37:42 0.0 (1 restrained)
## 6 -3738.361 437.305 1670 09:37:42 0.0
## 7 -3737.767 329.934 1670 09:37:42 0.0
## 8 -3737.397 405.785 1670 09:37:42 0.0
## 9 -3737.357 358.073 1670 09:37:42 0.0
## 10 -3737.312 390.537 1670 09:37:42 0.0
## 11 -3737.301 369.058 1670 09:37:42 0.0
## 12 -3737.293 383.625 1670 09:37:42 0.0
## 13 -3737.291 373.863 1670 09:37:42 0.0
## 14 -3737.289 380.468 1670 09:37:42 0.0
## 15 -3737.289 376.024 1670 09:37:42 0.0
## 16 -3737.288 379.026 1670 09:37:42 0.0
## 17 -3737.288 377.003 1670 09:37:42 0.0
## 18 -3737.288 378.369 1670 09:37:42 0.0asrmod(mod8)## $loglik
## [1] -3737.288
##
## $AIC
## [1] 7482.577
## attr(,"parameters")
## [1] 4
##
## $BIC
## [1] 7504.259
## attr(,"parameters")
## [1] 4summary(mod8)$varcomp## component std.error z.ratio bound %ch
## rep 3.828827e-05 7.483575e-06 5.116308 B 0
## row:col!R 3.783692e+02 7.395356e+01 5.116308 P 0
## row:col!row!cor 9.544919e-01 9.285171e-03 102.797448 U 0
## row:col!col!cor 2.858123e-01 2.444813e-02 11.690557 U 0wald(mod8)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 217183 574.00 <2e-16 ***
## cnitro 1 366 0.97 0.3256
## cnitro.sq 1 3 0.01 0.9249
## residual (MS) 378
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod8))
mod9 <- asreml(yield~ 1+cnitro+cnitro.sq+lin(row),
random= ~rep+spl(row),
residual=~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:43 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3630.728 26.7968 1669 09:37:44 0.0 (1 restrained)
## 2 -3532.309 25.1043 1669 09:37:44 0.0
## 3 -3465.017 25.4771 1669 09:37:44 0.0
## 4 -3444.899 27.0167 1669 09:37:44 0.0
## 5 -3441.556 28.2771 1669 09:37:44 0.0
## 6 -3441.455 28.5690 1669 09:37:44 0.0
## 7 -3441.452 28.6285 1669 09:37:44 0.0asrmod(mod9)## $loglik
## [1] -3441.452
##
## $AIC
## [1] 6892.904
## attr(,"parameters")
## [1] 5
##
## $BIC
## [1] 6920.004
## attr(,"parameters")
## [1] 5summary(mod9)$varcomp## component std.error z.ratio bound %ch
## rep 0.2028426 0.39285386 0.5163308 P 0.2
## spl(row) 9.3422237 2.71474268 3.4412925 P 0.0
## row:col!R 28.6285460 1.28083295 22.3515065 P 0.0
## row:col!row!cor 0.5088928 0.02007327 25.3517626 U 0.1
## row:col!col!cor 0.1488067 0.02590938 5.7433552 U 0.1wald(mod9)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1203272 42031 < 2.2e-16 ***
## cnitro 1 2532 88 < 2.2e-16 ***
## cnitro.sq 1 251 9 0.003052 **
## lin(row) 1 77066 2692 < 2.2e-16 ***
## residual (MS) 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod9))
plot(mod9)
mod10 <- asreml(yield~ 1+cnitro+cnitro.sq+lin(row),
random= ~rep+spl(row)+outl,
residual=~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:46 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3632.147 26.7959 1669 09:37:46 0.1 (2 restrained)
## 2 -3532.485 25.1011 1669 09:37:46 0.1 (1 restrained)
## 3 -3464.956 25.4798 1669 09:37:46 0.1 (1 restrained)
## 4 -3444.886 27.0197 1669 09:37:46 0.1 (1 restrained)
## 5 -3441.555 28.2779 1669 09:37:46 0.1 (1 restrained)
## 6 -3441.455 28.5691 1669 09:37:46 0.1
## 7 -3441.452 28.6286 1669 09:37:46 0.1mod10<- update(mod10)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:47 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3441.452 28.6396 1669 09:37:47 0.1
## 2 -3441.452 28.6403 1669 09:37:47 0.1asrmod(mod10)## $loglik
## [1] -3441.452
##
## $AIC
## [1] 6894.903
## attr(,"parameters")
## [1] 6
##
## $BIC
## [1] 6927.423
## attr(,"parameters")
## [1] 6summary(mod10)$varcomp## component std.error z.ratio bound %ch
## outl 2.898196e-06 1.297424e-07 22.3380859 B 0
## rep 2.028607e-01 3.927744e-01 0.5164815 P 0
## spl(row) 9.344648e+00 2.714414e+00 3.4426027 P 0
## row:col!R 2.864031e+01 1.282129e+00 22.3380859 P 0
## row:col!row!cor 5.089531e-01 2.007415e-02 25.3536564 U 0
## row:col!col!cor 1.487868e-01 2.590959e-02 5.7425402 U 0wald(mod10)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1203902 42035 < 2.2e-16 ***
## cnitro 1 2529 88 < 2.2e-16 ***
## cnitro.sq 1 251 9 0.003074 **
## lin(row) 1 77020 2689 < 2.2e-16 ***
## residual (MS) 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod10))
plot(mod10)
mod11 <- asreml(yield~ 1+cnitro+cnitro.sq+lin(row),
random= ~rep+spl(row)+outl+cnitro:spl(row)+cnitro.sq:spl(row),
residual=~ar1(row):ar1(col),
data=gridded_df_ordered,maxit=30)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:49 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -4106.350 23.9321 1669 09:37:50 0.3 (4 restrained)
## 2 -3891.952 22.8173 1669 09:37:50 0.3 (3 restrained)
## 3 -3713.606 23.2858 1669 09:37:50 0.3 (3 restrained)
## 4 -3601.118 24.6452 1669 09:37:51 0.3 (3 restrained)
## 5 -3522.435 25.4587 1669 09:37:51 0.3 (3 restrained)
## 6 -3466.290 26.0618 1669 09:37:51 0.3 (2 restrained)
## 7 -3442.216 26.8583 1669 09:37:52 0.3 (2 restrained)
## 8 -3433.955 27.1907 1669 09:37:52 0.3
## 9 -3433.943 27.2382 1669 09:37:52 0.3mod11<- update(mod11)## Model fitted using the gamma parameterization.
## ASReml 4.1.0 Tue Jul 21 09:37:53 2020
## Spline: design points closer than 0.0092 have been merged.
## LogLik Sigma2 DF wall cpu
## 1 -3433.943 27.2472 1669 09:37:54 0.3
## 2 -3433.943 27.2475 1669 09:37:54 0.3asrmod(mod11)## $loglik
## [1] -3433.943
##
## $AIC
## [1] 6883.887
## attr(,"parameters")
## [1] 8
##
## $BIC
## [1] 6927.247
## attr(,"parameters")
## [1] 8summary(mod11)$varcomp## component std.error z.ratio bound %ch
## outl 2.757257e-06 1.211081e-07 22.7669104 B 0
## rep 2.214044e-01 3.925147e-01 0.5640665 P 0
## spl(row) 9.282358e+00 2.702708e+00 3.4344659 P 0
## cnitro:spl(row) 1.102903e-08 4.844323e-10 22.7669104 B NA
## cnitro.sq:spl(row) 1.102903e-08 4.844323e-10 22.7669104 B NA
## row:col!R 2.724753e+01 1.196804e+00 22.7669104 P 0
## row:col!row!cor 4.819059e-01 2.088209e-02 23.0774740 U 0
## row:col!col!cor 1.611688e-01 2.597173e-02 6.2055493 U 0wald(mod11)## [0;34mWald tests for fixed effects.[0m
## [0;34mResponse: yield[0m
## [0;34mTerms added sequentially; adjusted for those above.[0m
##
## Df Sum of Sq Wald statistic Pr(Chisq)
## (Intercept) 1 1146164 42065 < 2.2e-16 ***
## cnitro 1 2736 100 < 2.2e-16 ***
## cnitro.sq 1 290 11 0.001114 **
## lin(row) 1 80351 2949 < 2.2e-16 ***
## residual (MS) 27
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(varioGram(mod11))
plot(mod11)
lrt(mod10,mod11)## Likelihood ratio test(s) assuming nested random models.
## (See Self & Liang, 1987)
##
## df LR-statistic Pr(Chisq)
## mod11/mod10 2 15.017 0.0001904 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit11<- fitted(mod11)
fit7 <- fitted(mod7)
plot(fit11,fit7)
abline(0,1)
cor(fit11,fit7)## [1] NAComparison
Contour plots of fitted values
Raw data
Intercept by GWR
By ASReml
By Bayesian model
Contour plots of predicted values
Predicted yield (quintals/ha) for an overall medium nitrogen treatment of 75.4 kg/ha.
Prediction by GWR
Prediction by Bayesian model