#' Model age-mixing pattern.
#'
#' Models the age-mixing pattern for the population simulated in Simpact. The
#' user specifies a time point, time window, and age group for which they would
#' like to obtain a summary of the age-mixing pattern. The user also specifies
#' whether only relationships that started within the time window should be
#' used.
#'
#' The \code{pattern.modeller} function uses the \code{link[nlme]{lme}}
#' function to build a linear mixed effects model regressing the partner's age
#' at the time the relationship started on his/her own age at the beginning of
#' the relationship. The models are stratified by gender. Each person can have
#' more than one relationship so a random intercept at the level of the person
#' is added to the model.
#'
#' The models also explicitly account for heteroskedastic variance, since the
#' variance in partner ages tends to grow for older ages.
#'
#' The function produces a list containing two elements. The first is a
#' subsetted dataframe of relationships that meet the inclusion criteria
#' specified in the function's arguments. This dataframe includes
#' population-average partner age predictions from the model. The second element
#' is a dataframe that contains key age-mixing pattern outputs from the model.
#' These outputs can be found in the second dataframe:
#'
#' \describe{
#' \item{slope}{Also known as the beta-coefficient, this is the
#' change in partner age for each year increase in age. It represents how fast
#' the average age differences grow in the population.}
#' \item{intercept}{This is the average partner age for the first age in the
#' specified interval.}
#' \item{power}{This represents how fast the variance in partner ages grows with
#' increasing age.}
#' \item{lowerpower}{This is the lower limit of the 95%
#' confidence interval for the power coefficient.}
#' \item{upperpower}{This is the upper limit of the 95% confidence interval for
#' the power coefficient.
#' }
#' \item{bvar}{This is the between-subject variance from the model.}
#' \item{wvar}{Also known as the residual error, this represents the amount of
#' within-subject variance from the model.}
#' }
#'
#' @param dataframe The dataframe that is produced by \code{\link{agemix.df.maker}}
#' @param agegroup Boundaries of the age group that should be retained, e.g.
#' c(15, 30). The interval is closed on the left and open on the right.
#' @param timepoint Point in time during the simulation to be used in the
#' calculation.
#' @param timewindow The length of time before the timepoint for which
#' relationships should be included, e.g. 1, representing one year before the
#' timepoint. This should be a whole number.
#' @param start This is a logical indicating that only relationships starting
#' after the beginning of the window should be used. If start = FALSE
#' relationships could start before the time window. This is the default.
#'
#' @return returns a list with two elements: 1. Dataframe of relationships
#' inputted into model, and 2. dataframe of model outputs.
#'
#' @examples
#' data(persreldf)
#' agemixpatdat <- pattern.modeller(dataframe = persreldf, agegroup = c(15, 30),
#' timewindow = 1, timepoint = 30)
#'
#' @importFrom magrittr %>%
#' @import dplyr
#' @export
#
#@importFrom nlme lme
#@importFrom nlme predict
#@importFrom nlme intervals
#@importFrom nlme VarCorr
pattern.modeller <- function(dataframe,
agegroup,
timepoint,
timewindow,
start = FALSE) {
#Warnings
if (!is.data.frame(dataframe)) {
stop("dataframe wrong type")
}
if (length(agegroup) != 2) {
stop("Need an upper and lower bound")
}
if (timepoint < 1) {
stop("Time must be at least 1 year")
}
if (timewindow < 0) {
stop("Window must be a whole number")
}
#Constants
time <- timepoint
window <- timepoint - timewindow
lwrage <- agegroup[1]
uprage <- agegroup[2]
if (start == TRUE) {
# This only includes relationships that started
# in the time window
men <- dataframe %>%
dplyr::mutate(age = time - TOB) %>%
dplyr::filter(episodeorder == 1 &
(FormTime <= time & FormTime >= window) &
age >= lwrage &
age < uprage &
Gender == "male" &
TOD > time) %>%
dplyr::mutate(agerelform0 = agerelform - lwrage)
women <- dataframe %>%
dplyr::mutate(age = time - TOB) %>%
dplyr::filter(episodeorder == 1 &
(FormTime <= time & FormTime >= window) &
age >= lwrage &
age < uprage &
Gender == "female" &
TOD > time) %>%
dplyr::mutate(agerelform0 = agerelform - lwrage)
} else {
# This includes all relationships that were ongoing
# at somepoint during the time window, but may have
# started long before the time window.
men <- dataframe %>%
dplyr::mutate(age = time - TOB) %>%
dplyr::filter(FormTime <= time &
DisTime > window &
age >= lwrage &
age < uprage &
Gender == "male" &
TOD > time) %>%
dplyr::distinct(ID, relid, .keep_all = TRUE) %>%
dplyr::mutate(agerelform0 = agerelform - lwrage)
women <- dataframe %>%
dplyr::mutate(age = time - TOB) %>%
dplyr::filter(FormTime <= time &
DisTime > window &
age >= lwrage &
age < uprage &
Gender == "female" &
TOD > time) %>%
dplyr::distinct(ID, relid, .keep_all = TRUE) %>%
dplyr::mutate(agerelform0 = agerelform - lwrage)
}
# malemodel <- lme(pagerelform ~ agerelform0,
# data = men,
# random = ~1 | ID,
# weights = varPower(value = 0.5, form = ~agerelform0 + 1),
# method = "REML")
#
# femalemodel <- lme(pagerelform ~ agerelform0,
# data = women,
# random = ~1 | ID,
# weights = varPower(value = 0.5, form = ~agerelform0 + 1),
# method = "REML")
#
# men$pred <- predict(malemodel, men, level = 0)
# women$pred <- predict(femalemodel, women, level = 0)
comb <- bind_rows(men, women)
modoutput <- matrix(nrow = 1, ncol = 14) %>%
as.data.frame()
colnames(modoutput) <- c("slopem", "slopew", "interceptm", "interceptw",
"powerm", "lowerpowerm", "upperpowerm", "powerw",
"lowerpowerw", "upperpowerw", "bvarm", "bvarw",
"wvarm", "wvarw")
#
# # Extract slope/Beta-coefficent of the fixed effects from model
# modoutput$slopem <- malemodel$coefficients$fixed[2]
# modoutput$slopew <- femalemodel$coefficients$fixed[2]
#
# # Extract population intercepts
# modoutput$interceptm <- malemodel$coefficients$fixed[1]
# modoutput$interceptw <- femalemodel$coefficients$fixed[1]
#
# # Extract power
# modoutput$powerm <- (attributes(malemodel$apVar)$Pars["varStruct.power"])
# #modoutput$lowerpowerm <- intervals(malemodel)$varStruct[,1]
# #modoutput$upperpowerm <- intervals(malemodel)$varStruct[,3]
#
# modoutput$powerw <- (attributes(femalemodel$apVar)$Pars["varStruct.power"])
# #modoutput$lowerpowerw <- intervals(femalemodel)$varStruct[, 1]
# #modoutput$upperpowerw <- intervals(femalemodel)$varStruct[, 3]
#
# # Extract between-participant variance
# modoutput$bvarm <- VarCorr(malemodel)[1] %>%
# as.numeric()
#
# modoutput$bvarw <- VarCorr(femalemodel)[1] %>%
# as.numeric()
#
# # Extract residual variance
# modoutput$wvarm <- VarCorr(malemodel)[2] %>%
# as.numeric()
#
# modoutput$wvarw <- VarCorr(femalemodel)[2] %>%
# as.numeric()
agemix.pieces <- list(comb, modoutput)
return(agemix.pieces)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.