## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, include = FALSE---------------------------------------------------
library(IndexWizard)

## -----------------------------------------------------------------------------
tn <- c("RZM", "RZN", "RZEo", "RZEn", "RZR", "RZKm", "RZKd", "RZH", "RZC", "RZS")

w_old <- c(0.45, 0.2, 0.15, 0, 0.1, 0.03, 0, 0, 0, 0.07)
names(w_old) <- tn; w_old

w_new <- c(0.36, 0.18, 0, 0.15, 0.07, 0.015, 0.015, 0.18, 0.03, 0)
names(w_new) <- tn; w_old


## -----------------------------------------------------------------------------
G <- matrix(
  c(1.0,0.13,0.13,0.07,-0.15,0.11,0.07,0.09,-0.02,0.04,
    0.13,1.0,0.23,0.28,0.43,0.25,0.22,0.78,0.13,0.46,
    0.13,0.23,1.0,0.92,0.02,0.09,-0.05,0.25,-0.1,0.19,
    0.07,0.28,0.92,1.0,0.06,0.08,-.03,0.31,-.1,0.25,
    -0.15,0.43,.02,0.06,1.0,0.32,0.19,0.41,0.04,0.15,
    0.11,0.25,0.09,0.08,0.32,1.0,0.0,0.25,0.04,0.13,
    0.07,0.22,-.05,-.03,0.19,0,1.0,0.23,0.05,0.10,
    0.09,0.78,0.25,0.31,0.41,0.25,0.23,1.0,0.1,0.57,
    -.02,0.13,-.10,-.1,0.04,0.04,0.05,0.1,1.0,0.02,
    0.04,0.46,0.19,0.25,0.15,0.13,0.10,0.57,0.02,1.0)
  ,byrow = TRUE, nrow = length(tn), ncol = length(tn), dimnames = list(tn, tn)
)
G <- G*12^2
G

## -----------------------------------------------------------------------------
r2 <- c(0.743, 0.673, 0.638, 0.717, 0.541, 0.635, 0.604, 0.720, 0.499, 0.764)
names(r2) <- tn; r2


## -----------------------------------------------------------------------------
r2_old <- r2[w_old != 0]
r2_new <- r2[w_new != 0]

## -----------------------------------------------------------------------------
deltautr <- c(0.28401392, 0.21637703, 0.17932963, 0.09986764, 0.08767239, 0.13273939)
names(deltautr) <- tn[w_old>0]


## -----------------------------------------------------------------------------
h2 <- c(0.314, 0.090, 0.194, 0.194, 0.013, 0.049, 0.033, 0.061, 0.014, 0.273)
names(h2) <- tn

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

H <- matrix (
  c(1.00,0.06,0.06,-.20,0.05,0.03,
    0.06,1.00,0.14,0.40,0.20,0.46,
    0.06,0.14,1.00,-.03,0.03,0.15,
    -.20,0.40,-.03,1.00,0.30,0.13,
    0.05,0.20,0.03,0.30,1.00,0.11,
    0.03,0.46,0.15,0.13,0.11,1.00),
  nrow=6, ncol=6,
  dimnames = list(tn[w_old != 0], tn[w_old != 0]))
H <- H * 144
H

## ----eval=FALSE, include=FALSE------------------------------------------------
#  library(IndexWizard)

## -----------------------------------------------------------------------------
res <- SelInd(
  w = w_old,
  G = G,
  r2 = r2_old
)


## -----------------------------------------------------------------------------
res <- SelInd(w = w_old,  G = G,  r2 = r2_old, verbose = FALSE)

## -----------------------------------------------------------------------------
summary(res)

## -----------------------------------------------------------------------------
res

## -----------------------------------------------------------------------------
res$w
res$b_scaled

## -----------------------------------------------------------------------------
res_old <- SelInd(w = w_old, G = G,  r2 = r2_old, verbose = FALSE)
res_new <- SelInd(w = w_new, G = G,  r2 = r2_new, verbose = FALSE)

## -----------------------------------------------------------------------------
round(res_old$b,2)
round(res_new$b,2)

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

ind <- c(RZM = 130, RZN = 110, RZEo = 105, RZEn = 106, RZR = 95,
         RZKm = 100, RZKd = 101, RZH = 115,  RZC = 90,  RZS = 120)

## -----------------------------------------------------------------------------
t(res_old$b_scaled) %*% ind[names(res_old$b_scaled)]

## -----------------------------------------------------------------------------
t(res_new$b_scaled) %*% ind[names(res_new$b_scaled)]

## -----------------------------------------------------------------------------
res_old <- SelInd(w = w_old, G = G,  r2 = r2_old, verbose = FALSE)
round(res_old$d_G_exp_scaled, 2)

## -----------------------------------------------------------------------------
h2
res_old <- SelInd(w = w_old, G = G,  r2 = r2_old, h2 = h2, verbose = FALSE)
round(res_old$d_P_exp_scaled, 2)

## -----------------------------------------------------------------------------
res_old <- SelInd(w = w_old, G = G,  r2 = r2_old, h2 = h2, i = 0.1, verbose = FALSE)
round(res_old$dG, 2)

## -----------------------------------------------------------------------------
round(res_old$d_G_exp, 2)
round(res_old$d_P_exp, 2)

## -----------------------------------------------------------------------------
res_old <- SelInd(w = w_old, G = G,  r2 = r2_old, h2 = h2, verbose = FALSE)
res_new <- SelInd(w = w_new, G = G,  r2 = r2_new, h2 = h2, verbose = FALSE)


## -----------------------------------------------------------------------------
round(res_new$d_G_exp_scaled - res_old$d_G_exp_scaled,2)

## -----------------------------------------------------------------------------
round(res_new$d_P_exp_scaled - res_old$d_P_exp_scaled,2)

## -----------------------------------------------------------------------------
res_old <- SelInd(w = w_old, G = G,  r2 = r2_old, verbose = FALSE)

## -----------------------------------------------------------------------------
round(res_old$r_IP, 2)

## -----------------------------------------------------------------------------
round(res_old$r_IH, 2)

## -----------------------------------------------------------------------------
round(res_old$del_d_scaled, 2)

## -----------------------------------------------------------------------------
res_old <- SelInd(w = w_old[w_old>0], G = G[w_old>0,w_old>0],  r2 = r2_old, verbose = FALSE)
round(deltautr - res_old$d_G_exp_scaled, 3)


## -----------------------------------------------------------------------------
res_old <- SelInd(
  w = w_old[w_old>0],
  G = G[w_old>0, w_old>0],
  r2 = r2_old,
  d_G_obs = deltautr,
  verbose = FALSE)
round(res_old$w_real - res_old$w, 2)
round(res_old$b_real - res_old$b_scaled, 2)

## -----------------------------------------------------------------------------
res_old_corRes <- SelInd(
  w = w_old[w_old>0],
  G = G[w_old>0, w_old>0],
  r2 = r2_old,
  d_G_obs = deltautr,
  H = H,
  verbose = FALSE)

## -----------------------------------------------------------------------------
round(res_old$w_real - res_old$w, 2) # uncorrelated residuals
round(res_old_corRes$w_real - res_old_corRes$w, 2) # correlated residuals