tests/testthat/test-latent.R

In Infusion: Inference Using Simulation

if (FALSE) {
  # this cannot work. Need to handle individual-levle vs pop-level latent vars distinctly
  cat(crayon::yellow("test reftable:\n"))
  
  rmidPO <- function(mu, sdG,sdE, sample.size=ssize) {
    transmitted <- rnorm(sd=sdG, n=sample.size) # offspring breeding value : latent variable
    OP <- rnorm(transmitted,sd=sdE, n=sample.size) # offspring phenotype
    midPP <- mu+rnorm(transmitted/2,sd=sdG/4+sdE/2, n=sample.size) # mid-parent phenotype 
    return(c(O=OP, midP=midPP, latent=transmitted))
  } # simulate means and variances of normal samples of size 'sample.size'
  set.seed(123)
  # pseudo-sample with stands for the actual data to be analyzed:  
  ssize <- 40
  (Sobs <- rmidPO(mu=1, sdG=0.5,sdE=1, sample.size=ssize) )
  # Uniform sampling in parameter space:
  npoints <- 600
  parsp <- data.frame(mu=runif(npoints,min=2.8,max=5.2),
                      s2=runif(npoints,min=0.4,max=2.4),sample.size=ssize)
  # Build simulation table:
  # set.seed(456) # this shows the considerable impact of the first reftable...
  simuls <- add_reftable(Simulate="myrnorm",par.grid=parsp)
  
  ## trivial projections which should produce an y=x regression:
  #mufit <- project("mu",stats=c("mean","var"),data=simuls, method="randomForest")
  #s2fit <- project("s2",stats=c("mean","var"),data=simuls, method="randomForest")
  mufit <- project("mu",stats=c("mean","var"),data=simuls)
  s2fit <- project("s2",stats=c("mean","var"),data=simuls)
  
  ## plots
  if (inherits(mufit,"HLfit")) mapMM(mufit,map.asp=1,
                                     plot.title=title(main="prediction of normal mean",xlab="mean",ylab="var"))
  if (inherits(s2fit,"HLfit")) mapMM(s2fit,map.asp=1,
                                     plot.title=title(main="prediction of normal var",xlab="mean",ylab="var"))
  
  ## apply projections on simulated statistics
  projectors <- list("MEAN"=mufit,"VAR"=s2fit)
  projectors <- list2env(projectors)
  corrSobs <- project(Sobs,projectors=projectors)
  corrSimuls <- project(simuls,projectors=projectors)
  
  
  # Infer surface:
  densv <- infer_SLik_joint(corrSimuls,stat.obs=corrSobs)
  # Usual workflow using inferred surface:
  slik_j <- MSL(densv) ## find the maximum of the log-likelihood surface
  #slik_j <- refine(slik_j,maxit=5, update_projectors=TRUE)
  slik_j <- refine(slik_j,maxit=5,update_projectors=TRUE)
  plot(slik_j)
  # etc:
  profile(slik_j,c(mu=4)) ## profile summary logL for given parameter value
  confint(slik_j,"mu") ## compute 1D confidence interval for given parameter
  plot1Dprof(slik_j,pars="s2",gridSteps=40) ## 1D profile
  
}
if (FALSE) { # example of distinct trainsample
  trainsample <- sample(nrow(slik_j$raw_data),1000)
  mufit <- project("mu",stats=c("mean","var"),data=slik_j$raw_data[trainsample,])
  s2fit <- project("s2",stats=c("mean","var"),data=slik_j$raw_data[trainsample,])
  
  ## plots
  if (inherits(mufit,"HLfit")) mapMM(mufit,map.asp=1,
                                     plot.title=title(main="prediction of normal mean",xlab="mean",ylab="var"))
  if (inherits(s2fit,"HLfit")) mapMM(s2fit,map.asp=1,
                                     plot.title=title(main="prediction of normal var",xlab="mean",ylab="var"))
  
  ## apply projections on simulated statistics
  projectors <- list("MEAN"=mufit,"VAR"=s2fit)
  projectors <- list2env(projectors)
  corrSobs <- project(Sobs,projectors=projectors)
  corrSimuls <- project(slik_j$raw_data[-trainsample,],projectors=projectors)
  
  
  # Infer surface:
  densv <- infer_SLik_joint(corrSimuls,stat.obs=corrSobs)
  # Usual workflow using inferred surface:
  slik_j <- MSL(densv) 
  plot(slik_j)
}


if (FALSE) { # example of reprojecting accumulated simulations
  remufit <- project("mu",stats=c("mean","var"),data=slik_j$raw_data, method="REML", train_cP_size=400, trainingsize=1000)
  res2fit <- project("s2",stats=c("mean","var"),data=slik_j$raw_data, method="REML", train_cP_size=400, trainingsize=1000)
  reprojectors <- list("MEAN"=remufit,"VAR"=res2fit)
  reprojectors <- list2env(reprojectors)
  recorrSobs <- project(Sobs,projectors=reprojectors)
  recorrSimuls <- project(slik_j$raw_data,projectors=reprojectors)
  # Infer surface:
  redensv <- infer_SLik_joint(recorrSimuls,stat.obs=recorrSobs)
  # Usual workflow using inferred suface:
  reslik_j <- MSL(redensv) ## find the maximum of the log-likelihood surface
  plot(reslik_j)
}

if (FALSE) { # 1-parameter example
  parsp1 <- data.frame(mu=4,
                      s2=runif(npoints,min=0.4,max=2.4),sample.size=ssize)
  # Build simulation table:
  simuls1 <- add_reftable(Simulate="myrnorm",par.grid=parsp1)
  s2fit1 <- project("s2",stats=c("mean","var"),data=simuls1)
  ## apply projections on simulated statistics
  projectors1 <- list("VAR"=s2fit1)
  projectors1 <- list2env(projectors1)
  corrSobs1 <- project(Sobs,projectors=projectors1)
  corrSimuls1 <- project(simuls1,projectors=projectors1)
  # Infer surface:
  densv1 <- infer_SLik_joint(corrSimuls1,stat.obs=corrSobs1)
  # Usual workflow using inferred surface:
  slik_j1 <- MSL(densv1) ## find the maximum of the log-likelihood surface
  slik_j1 <- refine(slik_j1,maxit=5, update_projectors=TRUE)
  plot(slik_j1)
  # etc:
  confint(slik_j1,"s2") ## compute 1D confidence interval for given parameter
  plot1Dprof(slik_j1,pars="s2",gridSteps=40) ## 1D profile
}