Last updated: 2019-11-21
Checks: 7 0
Knit directory: IITA_2019GS/
This reproducible R Markdown analysis was created with workflowr (version 1.5.0.9000). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20191121) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rproj.user/
Ignored: data/.DS_Store
Ignored: output/.DS_Store
Untracked files:
Untracked: data/GBSdataMasterList_31818.csv
Untracked: data/IITA_CleanedTrialData_72519.rds
Untracked: data/IITA_ExptDesignsDetected_72619.rds
Untracked: data/IITA_GBStoPhenoMaster_33018.csv
Untracked: data/NRCRI_GBStoPhenoMaster_40318.csv
Untracked: data/iita_blupsForCrossVal_72619.rds
Untracked: data/iita_blupsForCrossVal_outliersRemoved_73019.rds
Untracked: output/CrossVal_72719/
Untracked: output/CrossVal_73019/
Untracked: output/GEBV_BRNHT1_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_BRNHT1_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_DM_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_DM_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_DM_white_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_DM_white_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_DM_yellow_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_DM_yellow_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_HI_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_HI_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_IITA_OutliersRemovedFALSE_73119.csv
Untracked: output/GEBV_IITA_OutliersRemovedTRUE_73119.csv
Untracked: output/GEBV_MCMDS_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_MCMDS_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_PLTHT_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_PLTHT_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_RTSZ_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_RTSZ_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_TCHART_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_TCHART_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_logFYLD_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logFYLD_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_logFYLD_lowCMD_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logFYLD_nohav5_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logFYLD_propHAVpt5_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logFYLD_propHAVpt75_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logRTNO_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logRTNO_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_logRTNO_lowCMD_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logRTNO_nohav5_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logRTNO_propHAVpt5_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logRTNO_propHAVpt75_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logTOPYLD_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logTOPYLD_OutliersEqualTRUE_IITA_73119.rds
Untracked: output/GEBV_logTOPYLD_lowCMD_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logTOPYLD_nohav5_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logTOPYLD_propHAVpt5_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/GEBV_logTOPYLD_propHAVpt75_OutliersEqualFALSE_IITA_73119.rds
Untracked: output/PCscores_IITA_SamplesForGP_73119.rds
Untracked: output/sketch_of_crossVal_scheme.png
Untracked: workflowr_log.R
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | bfffb51 | wolfemd | 2019-11-21 | Publish the first set of analyses and files for IITA 2019 GS, |
Follow outlined GenomicPredictionChecklist and previous pipeline to process cassavabase data for ultimate genomic prediction.
This will cover IITA data, all years, all trials, downloaded from DB.
Using the Cassavabase search wizard:
Copy to data/DatabaseDownload_72419/
IITA’s entire DB download is pretty big. I used a remote machine cbsurobbins.biohpc.cornell.edu to do this processing quickly.
Note: GitHub filesize limit is 50 Mb, so this dataset cannot be shared there.
library(tidyverse); library(magrittr)
path<-"data/DatabaseDownload_72419/"
dbdata<-tibble(files=list.files(path = path)) %>%
mutate(Type=ifelse(grepl("metadata",files),"metadata","phenotype"),
files=map(files,~read.csv(paste0(path,.),
na.strings = c("#VALUE!",NA,".",""," ","-","\""),
stringsAsFactors = F) %>%
mutate_all(.,as.character)))
dbdata %<>%
filter(Type=="phenotype") %>%
select(-Type) %>%
unnest() %>%
left_join(dbdata %>%
filter(Type=="metadata") %>%
select(-Type) %>%
unnest() %>%
rename(programName=breedingProgramName,
programDescription=breedingProgramDescription,
programDbId=breedingProgramDbId) %>%
group_by(programName))
dim(dbdata)
dbdata %<>%
group_by(programName,locationName,studyYear,studyName,studyDesign,studyDescription,observationLevel) %>%
filter(observationLevel=="plot") %>%
nest(.key = TrialData)
dim(dbdata)521 observations 2272 trials total
WARNING: User input required! I create my own variable TrialType manually, using grepl and ifelse() expressions. By setting non-identifiable trials missing for TrialType, I can easily exclude them.
dbdata %<>%
mutate(TrialType=ifelse(grepl("CE|clonal|13NEXTgenC1",studyName,ignore.case = T),"CET",NA),
TrialType=ifelse(grepl("EC",studyName,ignore.case = T),"ExpCET",TrialType),
TrialType=ifelse(grepl("PYT",studyName,ignore.case = T),"PYT",TrialType),
TrialType=ifelse(grepl("AYT",studyName,ignore.case = T),"AYT",TrialType),
TrialType=ifelse(grepl("UYT",studyName,ignore.case = T),"UYT",TrialType),
TrialType=ifelse(grepl("geneticgain|gg|genetic gain",studyName,ignore.case = T),"GeneticGain",TrialType),
TrialType=ifelse(grepl("Cassava",studyName,ignore.case = T) & grepl("/",studyName),"GeneticGain",TrialType),
TrialType=ifelse((grepl("clonal evaluation trial",!grepl("genetic gain",studyDescription,ignore.case = T),
ignore.case = T)),"CET",TrialType),
TrialType=ifelse(grepl("preliminary yield trial",studyDescription,ignore.case = T),"PYT",TrialType),
TrialType=ifelse(grepl("Crossingblock|GS.C4.CB|cross",studyName) & is.na(TrialType),"CrossingBlock",TrialType),
TrialType=ifelse(grepl("NCRP",studyName) & is.na(TrialType),"NCRP",TrialType),
TrialType=ifelse(grepl("conservation",studyName) & is.na(TrialType),"Conservation",TrialType)) %>%
arrange(programName,studyYear,locationName) #%>% count(TrialType)Exclude non-identifiable trials
1948 trials from IITA.
Did this step on cbsurobbins, took lots of RAM
dbdata_long<-dbdata %>%
unnest() %>%
mutate(NOHAV=as.numeric(`plant.stands.harvested.counting.CO_334.0000010`)) %>%
select(-`plant.stands.harvested.counting.CO_334.0000010`) %>%
gather(Trait,Value,contains("CO_"),-NOHAV)
nrow(dbdata_long)/100000046.13M rows!
WARNING: User input required! Select the traits to be kept / analyzed, since the database download was indescriminant. In addition, manually give them abbreviations for convenience sake.
dbdata_long %<>%
select(Trait) %>%
distinct %>%
separate(Trait,c("TraitName","TraitCode"),".CO",remove = F,extra = 'merge') %>%
select(Trait,TraitName) %>%
distinct %>%
filter(grepl(paste0("cassava.mosaic.disease.severity.1.month|cassava.mosaic.disease.severity.3|",
"cassava.mosaic.disease.severity.6|cassava.mosaic.disease.severity.9|",
"dry.matter|total.carotenoid.by.chart.1.8|",
"plant.height.measurement.in.cm|first.apical.branch.height.measurement.in.cm|",
"fresh.shoot.weight.measurement.in.kg.per.plot|fresh.storage.root.weight.per.plot|",
"root.number.counting|storage.root.size.visual.rating.1.7"),
Trait,ignore.case = T)) %>%
filter(!grepl("Cassava.brown.streak.disease.leaf.severity.CO_334.0000036",Trait,ignore.case = T)) %>%
filter(!grepl("Cassava.brown.streak.disease.root.severity.CO_334.0000090",Trait,ignore.case = T)) %>%
filter(!grepl("marketable.root",Trait,ignore.case = T)) %>%
filter(!grepl("dry.matter.visual.rating.1.3",Trait,ignore.case = T)) %>%
mutate(TraitAbbrev=c("CMD1S","CMD3S","CMD6S","CMD9S",
"DMsg","DM",
"BRNHT1","SHTWT","RTWT","PLTHT","RTNO",
"RTSZ","TCHART")) %>%
inner_join(dbdata_long,.) %>%
rename(FullTraitName=Trait,
Trait=TraitAbbrev)
nrow(dbdata_long)/1000000Now only ~3.63M rows.
For each trait:
Deliberatiely leave out HI (calculate it manually after further QC)
dbdata_long %<>%
mutate(TraitType=ifelse(grepl("CBSD|CAD|CBB|CMD|CGM",Trait),"Disease",
ifelse(grepl("FYLD|RTWT|SHTWT|RTNO|DM|DMsg|RTSZ",Trait),"Yield","Misc")),
DiseaseScoreType=ifelse(TraitType=="Disease",
ifelse(grepl("S",Trait),"Severity","Incidence"),
NA))
dbdata_long %<>%
mutate(Value=as.numeric(Value),
Value=ifelse(TraitType=="Disease" & DiseaseScoreType=="Severity",
ifelse(Value<1 | Value>5,NA,Value),Value),
Value=ifelse(TraitType=="Disease" & DiseaseScoreType=="Incidence",
ifelse(Value<=0 | Value>1,NA,Value),Value),
Value=ifelse(Trait=="DM",
ifelse(Value>100 | Value<=0,NA,Value),Value),
Value=ifelse(Trait=="SPROUT",
ifelse(Value>1 | Value<=0,NA,Value),Value),
Value=ifelse(TraitType=="Yield",
ifelse(Value==0 | NOHAV==0 | is.na(NOHAV),NA,Value),Value),
NOHAV=ifelse(NOHAV==0,NA,NOHAV),
NOHAV=ifelse(NOHAV>42,NA,NOHAV),
Value=ifelse((Trait=="RTNO") & (!Value %in% 1:4000),NA,Value))Did this step on cbsurobbins, took lots of RAM
dbdata<-dbdata_long %>%
select(-FullTraitName,-TraitName,-TraitType,-DiseaseScoreType) %>%
spread(Trait,Value) %>%
mutate(DM=ifelse(is.na(DM) & !is.na(DMsg),DMsg,DM)) %>% # Fill in any missing DM scores with spec. grav-based scores
select(-DMsg)
rm(dbdata_long); gc()
nrow(dbdata)279595 obs left.
WARNING: User input required! At present, though cassavabase has the functionality to assign genotypes to phenotypes, the meta-information is in place; a breeding-program end task, perhaps. Instead, I rely on flat files, which I created over the years.
63K germplasmNames
library(tidyverse); library(magrittr)
gbs2phenoMaster<-dbdata %>%
select(germplasmName) %>%
distinct %>%
left_join(read.csv(paste0("data/",
"IITA_GBStoPhenoMaster_33018.csv"),
stringsAsFactors = F)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName) %>%
bind_rows(dbdata %>%
select(germplasmName) %>%
distinct %>%
left_join(read.csv(paste0("data/",
"NRCRI_GBStoPhenoMaster_40318.csv"),
stringsAsFactors = F)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName)) %>%
bind_rows(dbdata %>%
select(germplasmName) %>%
distinct %>%
left_join(read.csv("data/GBSdataMasterList_31818.csv",
stringsAsFactors = F) %>%
select(DNASample,FullSampleName) %>%
rename(germplasmName=DNASample)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName)) %>%
bind_rows(dbdata %>%
select(germplasmName) %>%
distinct %>%
mutate(germplasmSynonyms=ifelse(grepl("^UG",germplasmName,ignore.case = T),
gsub("UG","Ug",germplasmName),germplasmName)) %>%
left_join(read.csv("data/GBSdataMasterList_31818.csv",
stringsAsFactors = F) %>%
select(DNASample,FullSampleName) %>%
rename(germplasmSynonyms=DNASample)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName)) %>%
bind_rows(dbdata %>%
select(germplasmName) %>%
distinct %>%
mutate(germplasmSynonyms=ifelse(grepl("^TZ",germplasmName,
ignore.case = T),
gsub("TZ","",germplasmName),germplasmName)) %>%
left_join(read.csv("data/GBSdataMasterList_31818.csv",
stringsAsFactors = F) %>%
select(DNASample,FullSampleName) %>%
rename(germplasmSynonyms=DNASample)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName)) %>%
distinct %>%
left_join(read.csv("data/GBSdataMasterList_31818.csv",
stringsAsFactors = F) %>%
select(FullSampleName,OrigKeyFile,Institute) %>%
rename(OriginOfSample=Institute))
nrow(gbs2phenoMaster) #7866
gbs2phenoMaster %>% count(OriginOfSample)
# first, filter to just program-DNAorigin matches
germNamesWithGenos<-dbdata %>%
select(programName,germplasmName) %>%
distinct %>%
left_join(gbs2phenoMaster) %>%
filter(!is.na(FullSampleName))
nrow(germNamesWithGenos) # 7866
# program-germNames with locally sourced GBS samples
germNamesWithGenos_HasLocalSourcedGBS<-germNamesWithGenos %>% #count(OriginOfSample)
filter(programName==OriginOfSample) %>%
select(programName,germplasmName) %>%
semi_join(germNamesWithGenos,.) %>%
group_by(programName,germplasmName) %>% # select one DNA per germplasmName per program
slice(1) %>% ungroup()
nrow(germNamesWithGenos_HasLocalSourcedGBS) # 6816
# the rest (program-germNames) with GBS but coming from a different breeding program
germNamesWithGenos_NoLocalSourcedGBS<-germNamesWithGenos %>%
filter(programName==OriginOfSample) %>%
select(programName,germplasmName) %>%
anti_join(germNamesWithGenos,.) %>%
group_by(programName,germplasmName) %>% # select one DNA per germplasmName per program
slice(1) %>% ungroup()
nrow(germNamesWithGenos_NoLocalSourcedGBS) # 163
gbsForPhenos<-bind_rows(germNamesWithGenos_HasLocalSourcedGBS,
germNamesWithGenos_NoLocalSourcedGBS)
nrow(gbsForPhenos) # 6979
dbdata %<>%
left_join(gbsForPhenos) Compute harvest index after QC of RTWT and SHTWT above.
For calculating fresh root yield:
dbdata %<>%
mutate(PlotSpacing=ifelse(programName!="IITA",1,
ifelse(studyYear<2013,1,
ifelse(TrialType %in% c("CET","GeneticGain","ExpCET"),1,0.8))))
dbdata %<>%
group_by(programName,locationName,studyYear,studyName,studyDesign,studyDescription) %>%
summarize(MaxNOHAV=max(NOHAV, na.rm=T)) %>%
mutate(MaxNOHAV=ifelse(MaxNOHAV=="-Inf",NA,MaxNOHAV)) %>%
left_join(dbdata,.)WARNING: User input required! Only minor here. Depends on which disease traits are to be analyzed and which months-after-planting are recorded.
Compute season-wide mean (or if you wanted, AUDPC) after QC of trait values above.
WARNING: User input required! A few trials have variants on the most common / consensus locationName, so I have to fix them.
table(dbdata$locationName) # Showed some problem locationNames
dbdata %<>%
mutate(locationName=ifelse(locationName=="ibadan","Ibadan",locationName),
locationName=ifelse(locationName=="bwanga","Bwanga",locationName),
locationName=ifelse(locationName=="maruku","Maruku",locationName),
locationName=ifelse(locationName=="kasulu","Kasulu",locationName),
locationName=ifelse(locationName=="UKIRIGURU","Ukiriguru",locationName),
locationName=ifelse(grepl("NaCRRI",locationName),"Namulonge",locationName))
table(dbdata$locationName)WARNING: User input required! If I had preselected locations before downloading, this wouldn’t have been necessary.
dbdata %<>%
filter(locationName %in% c("Abuja","Ibadan","Ikenne","Ilorin","Jos","Kano",
"Malam Madori","Mokwa","Ubiaja","Umudike","Warri","Zaria"))
# count(TrialType,studyYear) %>% spread(studyYear,n) 238,673 x 65 obs remaining
Whatever design is reported to cassavabase cannot be universally trusted.
Examples: - Some trials appear to be complete blocked designs and the blockNumber is used instead of replicate, which is what most use. - Some complete block designs have nested, incomplete sub-blocks, others simply copy the “replicate” variable into the “blockNumber variable” - Some trials have only incomplete blocks but the incomplete block info might be in the replicate and/or the blockNumber column
One reason it might be important to get this right is that the variance among complete blocks might not be the same among incomplete blocks. If we treat a mixture of complete and incomplete blocks as part of the same random-effect (replicated-within-trial), we assume they have the same variance.
Also error variances might be heterogeneous among different trial-types (blocking scheme available) and/or plot sizes (maxNOHAV).
library(tidyverse);library(magrittr)
dbdata<-readRDS("data/IITA_CleanedTrialData_72519.rds") %>%
# custom selecting columns (not ideal)
select(programName,locationName,studyYear,trialType,TrialType,studyName,germplasmName,FullSampleName,
observationUnitDbId,replicate,blockNumber,
NOHAV,MaxNOHAV,
DM,RTNO,HI,FYLD,TOPYLD,MCMDS,
BRNHT1,PLTHT,TCHART,RTSZ) %>% #%$% summary(TCHART)
# extra QC of TCHART
mutate(TCHART=ifelse(TCHART %in% 1:8,TCHART,NA)) %>%
# custom create covariables for a custom trait requested by IYR
mutate(CMDcovar=MCMDS,
TCHARTcovar=TCHART) %>%
# custom variable selection for dplyr::gather()
gather(Trait,Value,DM:RTSZ) %>%
mutate(PropHAV=NOHAV/MaxNOHAV,
Value=ifelse(Trait %in% c("RTNO","FYLD","TOPYLD") & is.na(PropHAV),NA,Value)) %>%
# remove missing values
filter(!is.na(Value)) %>%
mutate(Value=ifelse(Trait %in% c("RTNO","FYLD","TOPYLD"),log(Value),Value),
Trait=ifelse(Trait %in% c("RTNO","FYLD","TOPYLD"),paste0("log",Trait),Trait)) %>%
# create explicitly nested experimental design variables
# intended for use in downstream analyses
mutate(yearInLoc=paste0(programName,"_",locationName,"_",studyYear),
trialInLocYr=paste0(yearInLoc,"_",studyName),
repInTrial=paste0(trialInLocYr,"_",replicate),
blockInRep=paste0(repInTrial,"_",blockNumber)) %>%
group_by(programName,locationName,studyYear,trialType,TrialType,studyName,Trait) %>%
nest(.key = TrialData) WARNING: User input required! In the code-chunk above, I do a’lot of customization. Columns are selected, an extra trait QC is added, and covariates for a custom-trait are created. Not ideal.
Code below is “standardized” but ad hoc.
# Define complete blocks
dbdata %>%
mutate(Nobs=map_dbl(TrialData,~nrow(.)),
MaxNOHAV=map_dbl(TrialData,~unique(.$MaxNOHAV)),
Nrep=map_dbl(TrialData,~length(unique(.$replicate))),
Nblock=map_dbl(TrialData,~length(unique(.$blockInRep))),
Nclone=map_dbl(TrialData,~length(unique(.$germplasmName))),
# median number of obs per clone
medObsPerClone=map_dbl(TrialData,
~count(.,germplasmName) %$% round(median(n),1)),
# median number of obs per replicate
medObsPerRep=map_dbl(TrialData,
~count(.,replicate) %$% round(median(n),1)),
# Define complete block effects based on the "replicate" variable
CompleteBlocks=ifelse(Nrep>1 & medObsPerClone==Nrep & Nobs!=Nrep,
TRUE,FALSE),
CompleteBlocks=ifelse(Nrep>1 & medObsPerClone!=Nrep &
medObsPerClone>1 & Nobs!=Nrep,
TRUE,CompleteBlocks)) -> x
# Additional trials with imperfect complete blocks
x %>%
# Some complete blocks may only be represented by the "blockNumber" column
mutate(medBlocksPerClone=map_dbl(TrialData,
~select(.,blockInRep,germplasmName) %>%
# median number of blockInRep per clone
distinct %>%
count(germplasmName) %$%
round(median(n))),
# If CompleteBlocks==FALSE (complete blocks not detected based on replicate)
# and if more than half the clones are represented in more than one block based on the blockInRep variable
# Copy the blockInRep values into the repInTrial column
# Recompute Nrep
# and declare CompleteBlocks==TRUE
TrialData=ifelse(medBlocksPerClone>1 & CompleteBlocks==FALSE,
map(TrialData,~mutate(.,repInTrial=blockInRep)),TrialData),
Nrep=map_dbl(TrialData,~length(unique(.$repInTrial))),
CompleteBlocks=ifelse(medBlocksPerClone>1 & CompleteBlocks==FALSE,
TRUE,CompleteBlocks)) -> y
# Define incomplete blocks
y %>%
mutate(repsEqualBlocks=map_lgl(TrialData,
~all(.$replicate==.$blockNumber)),
NrepEqualNblock=ifelse(Nrep==Nblock,TRUE,FALSE),
medObsPerBlockInRep=map_dbl(TrialData,
~count(.,blockInRep) %$% round(median(n),1))) -> z
z %<>% # Define complete blocked trials with nested sub-blocks
mutate(IncompleteBlocks=ifelse(CompleteBlocks==TRUE & Nobs!=Nblock & Nblock>1 & medObsPerBlockInRep>1 & NrepEqualNblock==FALSE,TRUE,FALSE))
table(z$IncompleteBlocks)
z %<>% # Define clearly unreplicated (CompleteBlocks==FALSE & Nrep==1) trials with nested sub-blocks
mutate(IncompleteBlocks=ifelse(CompleteBlocks==FALSE & Nobs!=Nblock & Nblock>1 & medObsPerBlockInRep>1 & Nrep==1,TRUE,IncompleteBlocks))
table(z$IncompleteBlocks)
z %<>% # Define additional trials with incomplete blocks (blockInRep) where CompleteBlocks==FALSE but Nrep>1 and Nrep==Block
mutate(IncompleteBlocks=ifelse(CompleteBlocks==FALSE & IncompleteBlocks==FALSE &
Nobs!=Nblock & Nblock>1 & Nobs!=Nrep &
medObsPerBlockInRep>1 & Nrep>1 & NrepEqualNblock==TRUE,TRUE,IncompleteBlocks))
z %<>% # Last few cases (2 trials actually) where Nrep>1 and Nblock>1 and Nrep!=Nblock but CompleteBlocks==FALSE
mutate(IncompleteBlocks=ifelse(CompleteBlocks==FALSE & IncompleteBlocks==FALSE &
Nobs!=Nblock & Nobs!=Nrep &
medObsPerBlockInRep>1 & Nrep>1,TRUE,IncompleteBlocks))
R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] workflowr_1.5.0.9000 Rcpp_1.0.3 rprojroot_1.3-2
[4] digest_0.6.22 later_1.0.0 R6_2.4.1
[7] backports_1.1.5 git2r_0.26.1 magrittr_1.5
[10] evaluate_0.14 stringi_1.4.3 rlang_0.4.1
[13] fs_1.3.1 promises_1.1.0 whisker_0.4
[16] rmarkdown_1.17 tools_3.6.1 stringr_1.4.0
[19] glue_1.3.1 httpuv_1.5.2 xfun_0.11
[22] yaml_2.2.0 compiler_3.6.1 htmltools_0.4.0
[25] knitr_1.26