## First, run the algorithm without including time as a covariate.
orderedPair <- pairData[pairData$infectionDiffY > 0, ]
## Create a variable called snpClose that will define probable links
# (<3 SNPs) and nonlinks (>12 SNPs) all pairs with between 2-12 SNPs
# will not be used to train.
orderedPair$snpClose <- ifelse(orderedPair$snpDist < 3, TRUE,
ifelse(orderedPair$snpDist > 12, FALSE, NA))
table(orderedPair$snpClose)
## Running the algorithm
# NOTE should run with nReps > 1.
resGen <- nbProbabilities(orderedPair = orderedPair,
indIDVar = "individualID",
pairIDVar = "pairID",
goldStdVar = "snpClose",
covariates = c("Z1", "Z2", "Z3", "Z4"),
label = "SNPs", l = 1,
n = 10, m = 1, nReps = 1)
## Merging the probabilities back with the pair-level data
nbResultsNoT <- merge(resGen[[1]], orderedPair, by = "pairID", all = TRUE)
## Estimating the serial interval
# Using hierarchical clustering with a 0.05 absolute difference cutoff
estimateSI(nbResultsNoT, indIDVar = "individualID",
timeDiffVar = "infectionDiffY", pVar = "pScaled",
clustMethod = "hc_absolute", cutoff = 0.05, initialPars = c(2, 2))
# \donttest{
# Using all pairs
estimateSI(nbResultsNoT, indIDVar = "individualID",
timeDiffVar = "infectionDiffY", pVar = "pScaled",
clustMethod = "none", initialPars = c(2, 2))
# # Using a shifted gamma distribution:
# # not allowing serial intervals of less than 3 months (0.25 years)
estimateSI(nbResultsNoT, indIDVar = "individualID",
timeDiffVar = "infectionDiffY", pVar = "pScaled",
clustMethod = "hc_absolute", cutoff = 0.05,
initialPars = c(2, 2), shift = 0.25)
# # Using multiple cutoffs
estimateSI(nbResultsNoT, indIDVar = "individualID",
timeDiffVar = "infectionDiffY", pVar = "pScaled",
clustMethod = "hc_absolute", cutoff = c(0.025, 0.05), initialPars = c(2, 2))
# }
## Adding confidence intervals
# NOTE should run with bootSamples > 2.
estimateSI(nbResultsNoT, indIDVar = "individualID",
timeDiffVar = "infectionDiffY", pVar = "pScaled",
clustMethod = "hc_absolute", cutoff = 0.05,
initialPars = c(2, 2), shift = 0.25, bootSamples = 2)
Run the code above in your browser using DataLab