HEF diallel analysis
2025-12-16
Description and data exploration
This is an analysis of the USFS HEF diallel test in support of the publication of Graham Ford’s thesis chapter. The purpose is to only estimate selection efficiency and determine appropriate age for selection based on tree height.
The test consisted of two test sites containing 78 fullsib families from 13 parents along with selfed seed, clonal material, and open-pollinated families. Each test consisted of four replications with line plots of 8 trees per family/treatment.
DBH and Height scatterplots by age
The scatterplots were produced after data cleaning where rows with NA values were dropped. Only fullsib families remained in the dataset.
Scatterplot of DBH vs Height with fitted linear regression
Fitting the models
Both univariate and multivariate models were fit using data from ages 7, 17, and 40. Age 30 data were excluded from analysis because it was only a subset of heights.
asr.ht7 <- asreml(fixed = ht7 ~ plant,
random = ~ female + and(male),
residual = ~id(units),
equate.levels = c("female","male"),
workspace = 128e6, pworkspace = 128e6,
data = dat1)## ASReml Version 4.2 18/12/2025 11:57:14
## LogLik Sigma2 DF wall
## 1 -5364.690 10.49361 3184 11:57:15
## 2 -5362.246 10.51473 3184 11:57:15
## 3 -5362.196 10.52051 3184 11:57:15
## 4 -5362.195 10.51969 3184 11:57:15asr.ht17 <- asreml(fixed = ht17 ~ plant,
random = ~ female + and(male),
residual = ~id(units),
equate.levels = c("female","male"),
workspace = 128e6, pworkspace = 128e6,
data = dat1)## ASReml Version 4.2 18/12/2025 11:57:15
## LogLik Sigma2 DF wall
## 1 -7171.052 32.59001 3185 11:57:15 ( 1 restrained)
## 2 -7167.307 32.81832 3185 11:57:15
## 3 -7165.703 32.73463 3185 11:57:15
## 4 -7165.431 32.70103 3185 11:57:15
## 5 -7165.417 32.69297 3185 11:57:15
## 6 -7165.417 32.69243 3185 11:57:15asr.ht40 <- asreml(fixed = ht40 ~ plant,
random = ~ female + and(male),
residual = ~id(units),
equate.levels = c("female","male"),
workspace = 128e6, pworkspace = 128e6,
data = dat1)## ASReml Version 4.2 18/12/2025 11:57:15
## LogLik Sigma2 DF wall
## 1 -5371.215 43.97858 2234 11:57:15
## 2 -5368.856 44.09913 2234 11:57:15
## 3 -5368.792 44.13935 2234 11:57:15
## 4 -5368.790 44.13268 2234 11:57:15Variance component summaries for univariate models
summary(asr.ht7)$varcomp## component std.error z.ratio bound %ch
## female 0.3436064 0.1497815 2.294051 P 0
## units!R 10.5196899 0.2641496 39.824744 P 0summary(asr.ht17)$varcomp## component std.error z.ratio bound %ch
## female 0.4942632 0.2322348 2.128291 P 0
## units!R 32.6924299 0.8207837 39.830751 P 0summary(asr.ht40)$varcomp## component std.error z.ratio bound %ch
## female 1.443171 0.6455429 2.235593 P 0.1
## units!R 44.132682 1.3240292 33.332106 P 0.0Multivariate model and variance component summaries
asr.mult.ht <- asreml(fixed = cbind(ht7, ht17, ht40) ~ trait + plant,
random = ~ corgh(trait):female + and(corgh(trait):male) +
trait:cross +
trait:plant:rep +
trait:plant:rep:col +
trait:plant:rep:row +
trait:plant:rep:plot,
residual = ~id(units):corgh(trait),
equate.levels = c("female","male"),
workspace = 128e6, pworkspace = 128e6,
data = dat1)## ASReml Version 4.2 18/12/2025 11:57:15
## LogLik Sigma2 DF wall
## 1 -17889.26 1.0 8605 11:57:16 ( 10 restrained)
## 2 -16118.68 1.0 8605 11:57:16 ( 4 restrained)
## 3 -15878.04 1.0 8605 11:57:16 ( 2 restrained)
## 4 -15826.37 1.0 8605 11:57:16 ( 2 restrained)
## 5 -15803.88 1.0 8605 11:57:16
## 6 -15798.57 1.0 8605 11:57:16
## 7 -15797.99 1.0 8605 11:57:17
## 8 -15797.95 1.0 8605 11:57:17
## 9 -15797.95 1.0 8605 11:57:17summary(asr.mult.ht)$varcomp## component std.error z.ratio bound
## trait:female!trait!ht17:!trait!ht7.cor 0.6910173 0.20190369 3.422510 U
## trait:female!trait!ht40:!trait!ht7.cor 0.5775225 0.22883625 2.523737 U
## trait:female!trait!ht40:!trait!ht17.cor 0.9299036 0.07477589 12.435874 U
## trait:female!trait_ht7 0.3249366 0.15657625 2.075261 P
## trait:female!trait_ht17 0.4197836 0.21024687 1.996622 P
## trait:female!trait_ht40 1.5537357 0.69813448 2.225554 P
## trait:plant:rep 4.7227279 1.53348697 3.079731 P
## trait:cross 0.2332794 0.08374576 2.785566 P
## trait:plant:rep:col 0.1601629 0.04581124 3.496149 P
## trait:plant:rep:row 0.1954529 0.11685153 1.672661 P
## trait:plant:rep:plot 1.2829371 0.16049400 7.993676 P
## units:trait!R 1.0000000 NA NA F
## units:trait!trait!ht17:!trait!ht7.cor 0.6324416 0.01154162 54.796616 U
## units:trait!trait!ht40:!trait!ht7.cor 0.5243011 0.01580756 33.167755 U
## units:trait!trait!ht40:!trait!ht17.cor 0.7304598 0.01077645 67.782943 U
## units:trait!trait_ht7 7.9037281 0.21800389 36.254987 P
## units:trait!trait_ht17 23.5938567 0.62657234 37.655439 P
## units:trait!trait_ht40 41.6197828 1.29680363 32.094129 P
## %ch
## trait:female!trait!ht17:!trait!ht7.cor 0.0
## trait:female!trait!ht40:!trait!ht7.cor 0.0
## trait:female!trait!ht40:!trait!ht17.cor 0.0
## trait:female!trait_ht7 0.0
## trait:female!trait_ht17 0.0
## trait:female!trait_ht40 0.0
## trait:plant:rep 0.0
## trait:cross 0.1
## trait:plant:rep:col 0.1
## trait:plant:rep:row 0.1
## trait:plant:rep:plot 0.0
## units:trait!R 0.0
## units:trait!trait!ht17:!trait!ht7.cor 0.0
## units:trait!trait!ht40:!trait!ht7.cor 0.0
## units:trait!trait!ht40:!trait!ht17.cor 0.0
## units:trait!trait_ht7 0.0
## units:trait!trait_ht17 0.0
## units:trait!trait_ht40 0.0Fitting a multivariate model with starting values obtained from the univariate models does not significantly affect variance components. Allowing cross variance to vary with age “idh(trait):cross” slightly improved correlations of additive genetic variance across ages
Multivariate model with starting values and variance component summaries
# univariate starting values
vF <- c(ht7 = 0.3436064, ht17 = 0.4942632, ht40 = 1.443171) # Var(female)
vR <- c(ht7 = 10.5196899, ht17 = 32.6924299, ht40 = 44.132682) # Var(residual)
# starting correlations
rF <- c(r17_7 = 0.60, r40_7 = 0.50, r40_17 = 0.90) # genetic (female)
rR <- c(r17_7 = 0.60, r40_7 = 0.50, r40_17 = 0.70) # residual
# other random terms
vCross <- c(ht7 = 0.01, ht17 = 0.30, ht40 = 1.30) # idh(trait):cross
vRep <- 4.70
vCol <- 0.16
vRow <- 0.19
vPlot <- 1.30
# assemble start.values in the exact ASReml parameter order for this model ---
start_vals <- c(
unname(vF[c("ht7","ht17","ht40")]),
unname(rF[c("r17_7","r40_7","r40_17")]),
unname(vCross[c("ht7","ht17","ht40")]),
vRep, vCol, vRow, vPlot,
unname(vR[c("ht7","ht17","ht40")]),
unname(rR[c("r17_7","r40_7","r40_17")]))
asr.mv.B <- asreml(
fixed = cbind(ht7, ht17, ht40) ~ trait + plant,
random = ~corgh(trait):female + and(corgh(trait):male) +
idh(trait):cross +
trait:plant:rep +
trait:plant:rep:col +
trait:plant:rep:row +
trait:plant:rep:plot,
residual = ~ id(units):corgh(trait),
equate.levels = c("female","male"),
start.values = start_vals,
data = dat1,
workspace = 128e6,
pworkspace = 128e6)## ASReml Version 4.2 18/12/2025 11:57:17
## LogLik Sigma2 DF wall
## 1 -17742.86 1.0 8605 11:57:17 ( 11 restrained)
## 2 -16117.97 1.0 8605 11:57:18 ( 7 restrained)
## 3 -15903.73 1.0 8605 11:57:18 ( 3 restrained)
## 4 -15836.16 1.0 8605 11:57:18 ( 2 restrained)
## 5 -15802.64 1.0 8605 11:57:18 ( 1 restrained)
## 6 -15791.56 1.0 8605 11:57:18
## 7 -15788.88 1.0 8605 11:57:18
## 8 -15788.36 1.0 8605 11:57:18
## 9 -15788.31 1.0 8605 11:57:19
## 10 -15788.31 1.0 8605 11:57:19summary(asr.mv.B)$varcomp## component std.error z.ratio bound
## trait:female!trait!ht17:!trait!ht7.cor 0.67614755 0.20022794 3.3768891 U
## trait:female!trait!ht40:!trait!ht7.cor 0.58345462 0.22989564 2.5379108 U
## trait:female!trait!ht40:!trait!ht17.cor 0.98255538 0.07399380 13.2788885 U
## trait:female!trait_ht7 0.34321619 0.15587999 2.2017976 P
## trait:female!trait_ht17 0.41363231 0.21110226 1.9593931 P
## trait:female!trait_ht40 1.41046870 0.68550695 2.0575557 P
## trait:plant:rep 4.73997016 1.53872555 3.0804520 P
## trait:cross!trait_ht7 0.02267712 0.06647618 0.3411316 P
## trait:cross!trait_ht17 0.31119862 0.16359605 1.9022380 P
## trait:cross!trait_ht40 1.37508511 0.45000097 3.0557381 P
## trait:plant:rep:col 0.16329716 0.04610308 3.5420009 P
## trait:plant:rep:row 0.18817587 0.11532960 1.6316355 P
## trait:plant:rep:plot 1.29536344 0.16040888 8.0753849 P
## units:trait!R 1.00000000 NA NA F
## units:trait!trait!ht17:!trait!ht7.cor 0.63198941 0.01154124 54.7592219 U
## units:trait!trait!ht40:!trait!ht7.cor 0.52317805 0.01582589 33.0583576 U
## units:trait!trait!ht40:!trait!ht17.cor 0.72985136 0.01081224 67.5023180 U
## units:trait!trait_ht7 7.90744991 0.21810783 36.2547738 P
## units:trait!trait_ht17 23.48362791 0.62484999 37.5828248 P
## units:trait!trait_ht40 40.98502087 1.28195203 31.9707915 P
## %ch
## trait:female!trait!ht17:!trait!ht7.cor 0
## trait:female!trait!ht40:!trait!ht7.cor 0
## trait:female!trait!ht40:!trait!ht17.cor 0
## trait:female!trait_ht7 0
## trait:female!trait_ht17 0
## trait:female!trait_ht40 0
## trait:plant:rep 0
## trait:cross!trait_ht7 0
## trait:cross!trait_ht17 0
## trait:cross!trait_ht40 0
## trait:plant:rep:col 0
## trait:plant:rep:row 0
## trait:plant:rep:plot 0
## units:trait!R 0
## units:trait!trait!ht17:!trait!ht7.cor 0
## units:trait!trait!ht40:!trait!ht7.cor 0
## units:trait!trait!ht40:!trait!ht17.cor 0
## units:trait!trait_ht7 0
## units:trait!trait_ht17 0
## units:trait!trait_ht40 0Family mean heritability
This is the family mean heritability estimated from the multivariate family model.
Note: The heritability is not being used but was calculated out of interest. The heritability estimates from the univariate model and multivariate animal model are consistent with previous estimations and Graham’s work.
## Estimate SE trait
## hi...1 0.6608675 0.11323804 Height7
## hi...2 0.5451715 0.12675693 Height17
## hi...3 0.7389297 0.08995777 Height40Selection efficiency estimates
Selection efficiency was estimated as the ratio of the expected correlated response in the target-age trait from selection at an earlier age to the direct response from selection at the target age.
sel_eff_tbl## # A tibble: 3 × 6
## sel_age target_age rA_to40 acc_sel acc_target selection_efficiency
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 7 40 0.578 0.887 0.916 0.559
## 2 17 40 0.930 0.897 0.916 0.911
## 3 40 40 1 0.916 0.916 1Although early selection at age 7 captures just over half of the age-40 genetic gain, delaying selection to age 17 recovers nearly all (~90%) of the potential response at age 40, indicating that most of the late-age genetic signal is already expressed by mid-rotation.
When expressed per year, early selection is substantially more efficient despite lower total response: selection at age 7 achieves ~56% of the age-40 response but delivers the highest gain per year (~0.08), selection at age 17 captures ~91% of the response with a moderate gain per year (~0.05), and direct selection at age 40 yields the full response but the lowest gain per year (~0.03) due to the long cycle length.
Note: I need help figuring out how I can estimate the selection efficiency at ages 10 and 12 to see what happens to the response and if it would be better to wait until either of those ages instead of selecting at ages 7 or 17.