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

## ----ttest_ex1, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # collect the data for the boys in the vector x and for the girls in the
# # vector y
# x<-sesamesim$postnumb[which(sesamesim$sex==1)]
# y<-sesamesim$postnumb[which(sesamesim$sex==2)]
# # execute student's t-test
# ttest <- t_test(x,y, var.equal = TRUE)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. The names of the means are x and y.
# results <- bain(ttest, "x = y; x > y; x < y")
# # display the results
# results
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex2, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # collect the data for the boys in the vector x and for the girs in the
# # vector y
# x<-sesamesim$postnumb[which(sesamesim$sex==1)]
# y<-sesamesim$postnumb[which(sesamesim$sex==2)]
# # execute student's t-test
# ttest <- t_test(x,y, var.equal = FALSE)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. The names of the means are x and y.
# results <- bain(ttest, "x = y; x > y; x < y")
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex3, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # compare the pre with the post measurements
# ttest <- t_test(sesamesim$prenumb,sesamesim$postnumb,paired = TRUE)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Use difference to refer to the difference between means
# results <- bain(ttest, "difference=0; difference>0; difference<0")
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex4, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data se
# library(bain)
# # compare post measurements with the reference value 30
# ttest <- t_test(sesamesim$postnumb)
# # Check name of estimate
# set.seed(100)
# # test hypotheses with bain versus the reference value 30. Use x to refer to the mean
# results <- bain(ttest, "x=30; x>30; x<30")
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex5, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # collect the data for the boys in the vector x and for the girs in the
# # vector y
# x<-sesamesim$postnumb[which(sesamesim$sex==1)]
# y<-sesamesim$postnumb[which(sesamesim$sex==2)]
# # execute student's t-test
# ttest <- t_test(x,y, var.equal = TRUE)
# # compute the pooled within standard deviation using the variance of x
# # (ttest$v[1]) and y (ttest$v[2])
# pwsd <- sqrt(((length(x) -1) * ttest$v[1] + (length(y)-1) * ttest$v[2])/
# ((length(x) -1) + (length(y) -1)))
# # print pwsd in order to be able to include it in the hypothesis. Its value
# # is 12.60
# print(pwsd)
# # set a seed value
# set.seed(100)
# # test hypotheses (the means of boy and girl differ less than .2 * pwsd =
# # 2.52 VERSUS the means differ more than .2 * pwsd = 2.52) with bain
# # note that, .2 is a value for Cohen's d reflecting a "small" effect, that
# # is, the means differ less or more than .2 pwsd.
# # use x and y to refer to the means.
# results <- bain(ttest, "x - y > -2.52 & x - y < 2.52")
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex6, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # make a factor of variable site
# sesamesim$site <- as.factor(sesamesim$site)
# # execute an analysis of variance using lm() which, due to the -1, returns
# # estimates of the means per group
# anov <- lm(postnumb~site-1,sesamesim)
# # take a look at the estimated means and their names
# coef(anov)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain
# results <- bain(anov, "site1=site2=site3=site4=site5; site2>site5>site1>
# site3>site4")
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex7, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data se
# library(bain)
# # make a factor of variable site
# sesamesim$site <- as.factor(sesamesim$site)
# # Center the covariate. If centered the coef() command below displays the
# # adjusted means. If not centered the intercepts are displayed.
# sesamesim$prenumb <- sesamesim$prenumb - mean(sesamesim$prenumb)
# # execute an analysis of covariance using lm() which, due to the -1, returns
# # estimates of the adjusted means per group
# ancov <- lm(postnumb~site+prenumb-1,sesamesim)
# # take a look at the estimated adjusted means, the regression coefficient
# # of the covariate and their names
# coef(ancov)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain
# results <- bain(ancov, "site1=site2=site3=site4=site5;
#                         site2>site5>site1>site3>site4")
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex8, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # execute a multiple regression using lm()
# regr <- lm(postnumb ~ age + peabody + prenumb,sesamesim)
# # take a look at the estimated regression coefficients and their names
# coef(regr)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Note that standardize = FALSE denotes that the
# # hypotheses are in terms of unstandardized regression coefficients
# results<-bain(regr, "age = 0 & peab=0 & pre=0 ; age > 0 & peab > 0 & pre > 0"
# , standardize = FALSE)
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)
# # Since it is only meaningful to compare regression coefficients if they are
# # measured on the same scale, bain can also evaluate standardized regression
# # coefficients (based on the seBeta function by Jeff Jones and Niels Waller):
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Note that standardize = TRUE denotes that the
# # hypotheses are in terms of standardized regression coefficients
# results<-bain(regr, "age = peab = pre ; pre > age > peab",standardize = TRUE)
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex9, eval=FALSE----------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # make a factor of variable site
# sesamesim$site <- as.factor(sesamesim$site)
# # execute an analysis of variance using lm() which, due to the -1, returns
# # estimates of the means per group
# anov <- lm(postnumb~site-1,sesamesim)
# # take a look at the estimated means and their names
# coef(anov)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain
# results1 <- bain(anov, "site1=site2=site3=site4=site5; site2>site5>site1>
# site3>site4", fraction = 1)
# # display the results
# print(results1)
# # obtain the descriptives table
# summary(results1, ci = 0.95)
# # execute a sensitivity analysis. See Hoijtink, Mulder,
# # van Lissa, and Gu (2019) for elaborations.
# set.seed(100)
# results2 <- bain(anov, "site1=site2=site3=site4=site5; site2>site5>site1>
# site3>site4",fraction = 2)
# print(results2)
# set.seed(100)
# results3 <- bain(anov, "site1=site2=site3=site4=site5; site2>site5>site1>
# site3>site4",fraction = 3)
# print(results3)

## ----ttest_ex20, eval=FALSE---------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# 
# # load mice (multiple imputation of missing data)
# # inspect the mice help file to obtain further information
# library(mice)
# 
# # create missing values in four variables from the sesamesim data set
# misdat <- cbind(sesamesim$prenumb,sesamesim$postnumb,sesamesim$funumb,sesamesim$peabody)
# colnames(misdat) <- c("prenumb","postnumb","funumb","peabody")
# misdat <- as.data.frame(misdat)
# 
# set.seed(1)
# pmis <- .80
# #
# for (i in 1:240){
#   uni<-runif(1)
#   if (pmis < uni) {
#     misdat$funumb[i]<-NA
#   }
#   uni<-runif(1)
#   if (pmis < uni) {
#     misdat$prenumb[i]<-NA
#     misdat$postnumb[i]<-NA
#   }
#   uni<-runif(1)
#   if (pmis < uni) {
#     misdat$peabody[i]<-NA
#   }
# }
# # print data summaries - note the missing values (NAs)
# summary(misdat)
# 
# # use mice to create 1000 imputed data matrices. Note that, the approach used
# # below # is only one manner in which mice can be instructed. Many other
# # options are available.
# M <- 1000
# out <- mice(data = misdat, m = M, seed=999, meth=c("norm","norm","norm","norm"), diagnostics = FALSE, printFlag = FALSE)
# 
# # create vectors in which 1000 fits and complexities can be stored for each
# # of two hypotheses and their complement
# fits1 <- vector("numeric",1000)
# compls1 <- vector("numeric",1000)
# fits2 <- vector("numeric",1000)
# compls2 <- vector("numeric",1000)
# fitscomplement <- vector("numeric",1000)
# complscomplement <- vector("numeric",1000)
# 
# # analyse each imputed data set with lm and bain, store the resulting
# # fits and complexities
# set.seed(100)
# for(i in 1:M) {
# regr <- lm(funumb~prenumb+postnumb,complete(out,i))
# result <- bain(regr,"prenumb=0 & postnumb=0;prenumb>0 & postnumb>0")
# fits1[i]<-result$fit$Fit[1]
# compls1[i]<-result$fit$Com[1]
# fits2[i]<-result$fit$Fit[2]
# compls2[i]<-result$fit$Com[2]
# fitscomplement[i]<-result$fit$Fit[4]
# complscomplement[i]<-result$fit$Com[4]
# }
# 
# # compute the Bayes factor from the imputed fits and complexities
# # see Equation 23 in Hoijtink, H., Gu, X., Mulder, J., and Rosseel, Y. (2019).
# # Computing Bayes Factors from Data with Missing Values.
# # Psychological Methods, 24, 253-268.
# 
# BF1u <- sum(fits1)/sum(compls1)
# BF2u <- sum(fits2)/sum(compls2)
# BFcomplementu <- sum(fitscomplement)/sum(complscomplement)
# 
# # compute PMPc
# PMPcH1 <- BF1u/(BF1u + BF2u + BFcomplementu)
# PMPcH2 <- BF2u/(BF1u + BF2u + BFcomplementu)
# PMPcHcomplement <- BFcomplementu/(BF1u + BF2u + BFcomplementu)

## ----ttest_ex10, eval=FALSE---------------------------------------------------
# # Load the bain and lavaan libraries. Visit www.lavaan.org for
# # lavaan mini-tutorials, examples, and elaborations
# library(bain)
# library(lavaan)
# 
# #  Specify and fit the confirmatory factor model
# model1 <- '
#     A =~ Ab + Al + Af + An + Ar + Ac
#     B =~ Bb + Bl + Bf + Bn + Br + Bc
# '
# # Use the lavaan sem function to execute the confirmatory factor analysis
# fit1 <- sem(model1, data = sesamesim, std.lv = TRUE)
# 
# # Inspect the parameter names
# coef(fit1)
# 
# # Formulate hypotheses, call bain, obtain summary stats
# hypotheses1 <-
# " A=~Ab > .6 & A=~Al > .6 & A=~Af > .6 & A=~An > .6 & A=~Ar > .6 & A=~Ac >.6 &
# B=~Bb > .6 & B=~Bl > .6 & B=~Bf > .6 & B=~Bn > .6 & B=~Br > .6 & B=~Bc >.6"
# set.seed(100)
# y <- bain(fit1,hypotheses1,fraction=1,standardize=TRUE)
# sy <- summary(y, ci = 0.95)
# 

## ----ttest_ex11, eval=FALSE---------------------------------------------------
# # Load the bain and lavaan libraries. Visit www.lavaan.org for
# # lavaan mini-tutorials, examples, and elaborations
# library(bain)
# library(lavaan)
# 
# # Specify and fit the latent regression model
# model2 <- '
#     A  =~ Ab + Al + Af + An + Ar + Ac
#     B =~ Bb + Bl + Bf + Bn + Br + Bc
# 
#     A ~ B + age + peabody
# '
# fit2 <- sem(model2, data = sesamesim, std.lv = TRUE)
# 
# # Inspect the parameter names
# coef(fit2)
# 
# # Formulate hypotheses, call bain, obtain summary stats
# 
# hypotheses2 <- "A~B > A~peabody = A~age = 0;
#                A~B > A~peabody > A~age = 0;
# A~B > A~peabody > A~age > 0"
# set.seed(100)
# y1 <- bain(fit2, hypotheses2, fraction = 1, standardize = TRUE)
# sy1 <- summary(y1, ci = 0.99)

## ----ttest_ex12, eval=FALSE---------------------------------------------------
# # Load the bain and lavaan libraries. Visit www.lavaan.org for
# # lavaan mini-tutorials, examples, and elaborations
# library(bain)
# library(lavaan)
# 
# # Specify A multiple group regression model
# model3 <- '
#     postnumb ~ prenumb + peabody
# '
# # Assign labels to the groups to be used when formulating
# # hypotheses
# sesamesim$sex <- factor(sesamesim$sex, labels = c("boy", "girl"))
# # Fit the multiple group regression model
# fit3 <- sem(model3, data = sesamesim, std.lv = TRUE, group = "sex")
# # Inspect the parameter names
# coef(fit3)
# 
# # Formulate and evaluate the hypotheses
# hypotheses3 <-
# "postnumb~prenumb.boy = postnumb~prenumb.girl & postnumb~peabody.boy = postnumb~peabody.girl;
#  postnumb~prenumb.boy < postnumb~prenumb.girl & postnumb~peabody.boy < postnumb~peabody.girl
# "
# set.seed(100)
# y1 <- bain(fit3, hypotheses3, standardize = TRUE)
# sy1 <- summary(y1, ci = 0.95)

## ----ttest_ex13, eval=FALSE---------------------------------------------------
# # load the bain package which includes the simulated sesamesim data se
# library(bain)
# # make a factor of variable site
# sesamesim$site <- as.factor(sesamesim$site)
# # execute an analysis of variance using lm() which, due to the -1, returns
# # estimates of the means per group
# anov <- lm(postnumb~site-1,sesamesim)
# # take a look at the estimated means and their names
# coef(anov)
# # collect the estimates means in a vector
# estimate <- coef(anov)
# # give names to the estimates in anov
# names(estimate) <- c("site1", "site2", "site3","site4","site5")
# # create a vector containing the sample size of each group
# # (in case of missing values in the variables used, the command
# # below has to be modified such that only the cases without
# # missing values are counted)
# ngroup <- table(sesamesim$site)
# # compute for each group the covariance matrix of the parameters
# # of that group and collect these in a list
# # for the ANOVA this is simply a list containing for each group the variance
# # of the mean note that, the within group variance as estimated using lm is
# # used to compute the variance of each of the means! See, Hoijtink, Gu, and
# # Mulder (2019) for further elaborations.
# var <- summary(anov)$sigma**2
# cov1 <- matrix(var/ngroup[1], nrow=1, ncol=1)
# cov2 <- matrix(var/ngroup[2], nrow=1, ncol=1)
# cov3 <- matrix(var/ngroup[3], nrow=1, ncol=1)
# cov4 <- matrix(var/ngroup[4], nrow=1, ncol=1)
# cov5 <- matrix(var/ngroup[5], nrow=1, ncol=1)
# covlist <- list(cov1, cov2, cov3, cov4,cov5)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Note that there are multiple groups
# # characterized by one mean, therefore group_parameters=0. Note that are no
# # joint parameters, therefore, joint_parameters=0.
# results <- bain(estimate,
# "site1=site2=site3=site4=site5; site2>site5>site1>site3>site4",
# n=ngroup,Sigma=covlist,group_parameters=1,joint_parameters = 0)
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex21, eval=FALSE---------------------------------------------------
# # Load the bain library
# library(bain)
# 
# # make a factor of variable site
# sesamesim$site <- as.factor(sesamesim$site)
# 
# # collect the estimated means in a vector
# estimates <- aggregate(sesamesim$postnumb,list(sesamesim$site),mean)[,2]
# # give names to the estimates
# names(estimates) <- c("site1", "site2", "site3","site4","site5")
# 
# # create a vector containing the sample sizes of each group
# ngroup <- table(sesamesim$site)
# 
# # compute the variances in each group and subquently the squared standard errors
# vars <- aggregate(sesamesim$postnumb,list(sesamesim$site),var)
# vargroup <- vars[,2]/ngroup
# 
# # collect the variances of the means in a covariance listcov1 <- matrix(vargroup[1], nrow=1, ncol=1)
# cov2 <- matrix(vargroup[2], nrow=1, ncol=1)
# cov3 <- matrix(vargroup[3], nrow=1, ncol=1)
# cov4 <- matrix(vargroup[4], nrow=1, ncol=1)
# cov5 <- matrix(vargroup[5], nrow=1, ncol=1)
# covlist <- list(cov1, cov2, cov3, cov4,cov5)
# 
# # set a seed value
# set.seed(100)
# # test hypotheses using bain. Note that there are multiple groups
# # characterized by one mean, therefore group_parameters=1. Note that
# # there are no joint parameters, therefore, joint_parameters=0.
# results <- bain(estimates,
#                 "site1=site2=site3=site4=site5;site2>site5>site1>site4>site3",
#                 n=ngroup,Sigma=covlist,group_parameters=1,joint_parameters = 0)
# print(results)
# 
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex19, eval=FALSE---------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # load the WRS2 package which renders the trimmed sample mean and
# # corresponding standard error
# library(WRS2)
# # make a factor of variable site
# sesamesim$site <- as.factor(sesamesim$site)
# # create a vector containing the sample sizes of each group
# # (in case of missing values in the variables used, the command
# # below has to be modified such that only the cases without
# # missing values are counted)
# ngroup <- table(sesamesim$site)
# # Compute the 20\% sample trimmed mean for each site
# estimates <- c(mean(sesamesim$postnumb[sesamesim$site == 1], tr = 0.2),
#                mean(sesamesim$postnumb[sesamesim$site == 2], tr = 0.2),
#                mean(sesamesim$postnumb[sesamesim$site == 3], tr = 0.2),
#                mean(sesamesim$postnumb[sesamesim$site == 4], tr = 0.2),
#                mean(sesamesim$postnumb[sesamesim$site == 5], tr = 0.2))
# # give names to the estimates
# names(estimates) <- c("s1", "s2", "s3","s4","s5")
# # display the estimates and their names
# print(estimates)
# # Compute the sample trimmed mean standard error for each site
# se <- c(trimse(sesamesim$postnumb[sesamesim$site == 1]),
#         trimse(sesamesim$postnumb[sesamesim$site == 2]),
#         trimse(sesamesim$postnumb[sesamesim$site == 3]),
#         trimse(sesamesim$postnumb[sesamesim$site == 4]),
#         trimse(sesamesim$postnumb[sesamesim$site == 5]))
# # Square the standard errors to obtain the variances of the sample
# # trimmed means
# var <- se^2
# # Store the variances in a list of matrices
# covlist <- list(matrix(var[1]),matrix(var[2]),
# matrix(var[3]),matrix(var[4]), matrix(var[5]))
# # set a seed value
# set.seed(100)
# # test hypotheses with bain
# results <- bain(estimates,"s1=s2=s3=s4=s5;s2>s5>s1>s3>s4",
# n=ngroup,Sigma=covlist,group_parameters=1,joint_parameters= 0)
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex14, eval=FALSE---------------------------------------------------
# # load the bain package which includes the simulated sesamesim data se
# library(bain)
# # make a factor of variable site
# sesamesim$site <- as.factor(sesamesim$site)
# # center the covariate. If centered the coef() command below displays the
# # adjusted means. If not centered the intercepts are displayed.
# sesamesim$prenumb <- sesamesim$prenumb - mean(sesamesim$prenumb)
# # execute an analysis of covariance using lm() which, due to the -1, returns
# # estimates of the adjusted means per group
# ancov2 <- lm(postnumb~site+prenumb-1,sesamesim)
# # take a look at the estimated adjusted means and their names
# coef(ancov2)
# # collect the estimates of the adjusted means and regression coefficient of
# # the covariate in a vector (the vector has to contain first the
# # adjusted means and next the regression coefficients of the covariates)
# estimates <- coef(ancov2)
# # assign names to the estimates
# names(estimates)<- c("v.1", "v.2", "v.3", "v.4","v.5", "pre")
# # compute the sample size per group
# # (in case of missing values in the variables used, the command
# # below has to be modified such that only the cases without
# # missing values are counted)
# ngroup <- table(sesamesim$site)
# # compute for each group the covariance matrix of the parameters of that
# # group and collect these in a list note that, the residual variance as
# # estimated using lm is used to compute these covariance matrices
# var <- (summary(ancov2)$sigma)**2
# # below, for each group, the covariance matrix of the adjusted mean and
# # covariate is computed
# # see Hoijtink, Gu, and Mulder (2019) for further explanation and elaboration
# cat1 <- subset(cbind(sesamesim$site,sesamesim$prenumb), sesamesim$site == 1)
# cat1[,1] <- 1
# cat1 <- as.matrix(cat1)
# cov1 <- var * solve(t(cat1) %*% cat1)
# #
# cat2 <- subset(cbind(sesamesim$site,sesamesim$prenumb), sesamesim$site == 2)
# cat2[,1] <- 1
# cat2 <- as.matrix(cat2)
# cov2 <- var * solve(t(cat2) %*% cat2)
# #
# cat3 <- subset(cbind(sesamesim$site,sesamesim$prenumb), sesamesim$site == 3)
# cat3[,1] <- 1
# cat3 <- as.matrix(cat3)
# cov3 <- var * solve(t(cat3) %*% cat3)
# #
# cat4 <- subset(cbind(sesamesim$site,sesamesim$prenumb), sesamesim$site == 4)
# cat4[,1] <- 1
# cat4 <- as.matrix(cat4)
# cov4 <- var * solve(t(cat4) %*% cat4)
# #
# cat5 <- subset(cbind(sesamesim$site,sesamesim$prenumb), sesamesim$site == 5)
# cat5[,1] <- 1
# cat5 <- as.matrix(cat5)
# cov5 <- var * solve(t(cat5) %*% cat5)
# #
# covariances <- list(cov1, cov2, cov3, cov4,cov5)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Note that there are multiple groups
# # characterized by one adjusted mean, therefore group_parameters=1 Note that
# # there is one covariate, therefore, joint_parameters=1.
# results2<-bain(estimates,"v.1=v.2=v.3=v.4=v.5;v.2 > v.5 > v.1 > v.3 >v.4;",
# n=ngroup,Sigma=covariances,group_parameters=1,joint_parameters = 1)
# # display the results
# print(results2)
# # obtain the descriptives table
# summary(results2, ci = 0.95)

## ----ttest_ex15, eval=FALSE---------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # estimate the means of three repeated measures of number knowledge
# # the 1 denotes that these means have to be estimated
# within <- lm(cbind(prenumb,postnumb,funumb)~1, data=sesamesim)
# # take a look at the estimated means and their names
# coef(within)
# # note that the model specified in lm has three dependent variables.
# # Consequently, the estimates rendered by lm are collected in a "matrix".
# # Since bain needs a named vector containing the estimated means, the [1:3]
# # code is used to select the three means from a matrix and store them in a
# # vector.
# estimate <- coef(within)[1:3]
# # give names to the estimates in anov
# names(estimate) <- c("pre", "post", "fu")
# # compute the sample size
# # (in case of missing values in the variables used, the command
# # below has to be modified such that only the cases without
# # missing values are counted)
# ngroup <- nrow(sesamesim)
# # compute the covariance matrix of the three means
# covmatr <- list(vcov(within))
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Note that there is one group, with three
# # means therefore group_parameters=3.
# results <- bain(estimate,"pre = post = fu; pre < post < fu", n=ngroup,
# Sigma=covmatr, group_parameters=3, joint_parameters = 0)
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex16, eval=FALSE---------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # make a factor of the variable sex
# sesamesim$sex <- factor(sesamesim$sex)
# # estimate the means of prenumb, postnumb, and funumb for boys and girls
# # the -1 denotes that the means have to estimated
# bw <- lm(cbind(prenumb, postnumb, funumb)~sex-1, data=sesamesim)
# # take a look at the estimated means and their names
# coef(bw)
# # collect the estimated means in a vector
# est1 <-coef(bw)[1,1:3] # the three means for sex = 1
# est2 <-coef(bw)[2,1:3] # the three means for sex = 2
# estimate <- c(est1,est2)
# # give names to the estimates in anov
# names(estimate) <- c("pre1", "post1", "fu1","pre2", "post2", "fu2")
# # determine the sample size per group
# # (in case of missing values in the variables used, the command
# # below has to be modified such that only the cases without
# # missing values are counted)
# ngroup<-table(sesamesim$sex)
# # cov1 has to contain the covariance matrix of the three means in group 1.
# # cov2 has to contain the covariance matrix in group 2
# # typing vcov(bw) in the console pane highlights the structure of
# # the covariance matrix of all 3+3=6 means
# # it has to be dissected in to cov1 and cov2
# cov1 <- c(vcov(bw)[1,1],vcov(bw)[1,3],vcov(bw)[1,5],vcov(bw)[3,1],
# vcov(bw)[3,3],vcov(bw)[3,5],vcov(bw)[5,1],vcov(bw)[5,3],vcov(bw)[5,5])
# cov1 <- matrix(cov1,3,3)
# cov2 <- c(vcov(bw)[2,2],vcov(bw)[2,4],vcov(bw)[2,6],vcov(bw)[4,2],
# vcov(bw)[4,4],vcov(bw)[4,6],vcov(bw)[6,2],vcov(bw)[6,4],vcov(bw)[6,6])
# cov2 <- matrix(cov2,3,3)
# covariance<-list(cov1,cov2)
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Note that there are multiple groups
# # characterized by three means, therefore group_parameters=3. Note that there
# # are no additional parameters, therefore, joint_parameters=0.
# results <-bain(estimate, "pre1 - pre2 = post1 - post2 = fu1 -fu2;
# pre1 - pre2 > post1 - post2 > fu1 -fu2"  , n=ngroup, Sigma=covariance,
# group_parameters=3, joint_parameters = 0)
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex17, eval=FALSE---------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # regression coefficients can only be mutually compared if the
# # corresponding predictors are measured on the same scale, therefore the
# # predictors are standardized
# sesamesim$age <- (sesamesim$age - mean(sesamesim$age))/sd(sesamesim$age)
# sesamesim$peabody <- (sesamesim$peabody - mean(sesamesim$peabody))/
# sd(sesamesim$peabody)
# sesamesim$prenumb <- (sesamesim$prenumb - mean(sesamesim$prenumb))/
# sd(sesamesim$prenumb)
# # estimate the logistic regression coefficients
# logreg <- glm(viewenc ~ age + peabody + prenumb, family=binomial,
# data=sesamesim)
# # take a look at the estimates and their names
# coef(logreg)
# # collect the estimated intercept and regression coefficients in a vector
# estimate <- coef(logreg)
# # give names to the estimates
# names(estimate) <- c("int", "age", "peab" ,"pre" )
# # compute the sample size. Note that, this is an analysis with ONE group
# # (in case of missing values in the variables used, the command
# # below has to be modified such that only the cases without
# # missing values are counted)
# ngroup <- nrow(sesamesim)
# # compute the covariance matrix of the intercept and regression coefficients
# covmatr <- list(vcov(logreg))
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Note that there is one group, with four
# # parameters (intercept plus three regression coefficients) therefore
# # group_parameters=4.
# results <- bain(estimate, "age = peab = pre; age > pre > peab", n=ngroup,
# Sigma=covmatr, group_parameters=4, joint_parameters = 0)
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----ttest_ex18, eval=FALSE---------------------------------------------------
# # load the bain package which includes the simulated sesamesim data set
# library(bain)
# # load the numDeriv package which will be used to compute the covariance
# # matrix of the adjusted group coefficients and regression coefficient of age
# # for the boys and the girls # using the estimates obtained using the data for
# # the boys AND the girls.
# library(numDeriv)
# # make a factor of the variable sex
# sesamesim$sex <- factor(sesamesim$sex)
# # center the covariate age
# sesamesim$age <- sesamesim$age - mean(sesamesim$age)
# # determine sample size per sex group
# # (in case of missing values in the variables used, the command
# # below has to be modified such that only the cases without
# # missing values are counted)
# ngroup <- table(sesamesim$sex)
# # execute the logistic regression, -1 ensures that the coefficients
# # for boys and girl are estimated adjusted for the covariate age
# anal <- glm(viewenc ~ sex + age - 1, family=binomial, data=sesamesim)
# # take a look at the estimates and their names
# coef(anal)
# # collect the estimates obtained using the data of the boys AND the
# # girls in a vector. This vector contains first the group specific
# # parameters followed by the regression coefficient of the covariate.
# estimates <-coef(anal)
# # give names to the estimates
# names(estimates) <- c("boys", "girls", "age")
# # use numDeriv to compute the Hessian matrix and subsequently the
# # covariance matrix for each of the two (boys and girls) groups.
# # The vector f should contain the regression coefficient of the group
# # at hand and the regression coefficient of the covariate.
# #
# # the first group
# data <- subset(cbind(sesamesim$sex,sesamesim$age,sesamesim$viewenc),
# sesamesim$sex==1)
# f <- 1
# f[1] <- estimates[1] # the regression coefficient of boys
# f[2] <- estimates[3] # the regression coefficient of age
# #
# # within the for loop below the log likelihood of the logistic
# # regression model is computed using the data for the boys
# logist1 <- function(x){
#   out <- 0
#   for (i in 1:ngroup[1]){
#     out <- out + data[i,3]*(x[1] + x[2]*data[i,2]) - log (1 +
#     exp(x[1] + x[2]*data[i,2]))
#   }
#   return(out)
# }
# hes1 <- hessian(func=logist1, x=f)
# # multiply with -1 and invert to obtain the covariance matrix for the
# # first group
# cov1 <- -1 * solve(hes1)
# #
# # the second group
# data <- subset(cbind(sesamesim$sex,sesamesim$age,sesamesim$viewenc),
# sesamesim$sex==2)
# f[1] <- estimates[2] # the regression coefficient of girls
# f[2] <- estimates[3] # the regression coefficient of age
# 
# # within the for loop below the log likelihood of the logistic
# # regression model is computed using the data for the girls
# logist2 <- function(x){
#   out <- 0
#   for (i in 1:ngroup[2]){
#     out <- out + data[i,3]*(x[1] + x[2]*data[i,2]) - log (1 +
#     exp(x[1] + x[2]*data[i,2]))
#   }
#   return(out)
# }
# hes2 <- hessian(func=logist2, x=f)
# # multiply with -1 and invert to obtain the covariance matrix
# cov2 <- -1 * solve(hes2)
# #
# #make a list of covariance matrices
# covariance<-list(cov1,cov2)
# #
# # set a seed value
# set.seed(100)
# # test hypotheses with bain. Note that there are multiple groups
# # characterized by one adjusted group coefficient,
# # therefore group_parameters=1. Note that there is one covariate,
# # therefore, joint_parameters=1.
# results <- bain(estimates, "boys < girls & age > 0", n=ngroup,
# Sigma=covariance, group_parameters=1, joint_parameters = 1)
# # display the results
# print(results)
# # obtain the descriptives table
# summary(results, ci = 0.95)

## ----eval = FALSE-------------------------------------------------------------
# library(bain)
# library(lavaan)
# model1 <- 'A =~ Ab + Al + Af + An + Ar + Ac
#            B =~ Bb + Bl + Bf + Bn + Br + Bc'
# fit1 <- sem(model1, data = sesamesim, std.lv = TRUE)

## ----eval = FALSE-------------------------------------------------------------
# hypotheses <- "(Ab, Al, Af, An, Ar, Ac) >.6 &
#                (Bb, Bl, Bf, Bn, Br, Bc) >.6"

## ----eval = FALSE-------------------------------------------------------------
# # Extract standardized parameter estimates (argument x)
# PE1 <- parameterEstimates(fit1, standardized = TRUE)
# # Identify which parameter estimates are factor loadings
# loadings <- which(PE1$op == "=~")
# # Collect the standardized "std.all" factor loadings "=~"
# estimates <- PE1$std.all[loadings]
# # Assign names to the standardized factor loadings
# names(estimates) <- c("Ab", "Al", "Af", "An", "Ar", "Ac",
#                       "Bb", "Bl", "Bf", "Bn", "Br", "Bc")
# 
# # Extract sample size (argument n)
# n <- nobs(fit1)
# 
# # Extract the standardized model parameter covariance matrix,
# # selecting only the rows and columns for the loadings
# covariance <- lavInspect(fit1, "vcov.std.all")[loadings, loadings]
# 
# # Give the covariance matrix the same names as the estimate
# dimnames(covariance) <- list(names(estimates), names(estimates))
# 
# # Put the covariance matrix in a list
# covariance <- list(covariance)
# 
# # Specify number of group parameters; this simply means that there
# # are 12 parameters in the estimate.
# group_parameters <- 12
# 
# # Then, the number of joint parameters. This is always 0 for evaluations of
# # SEM hypotheses:
# joint_parameters <- 0

## ----eval = FALSE-------------------------------------------------------------
# res <- bain(x = estimates,
#             hypothesis = hypotheses,
#             n = n,
#             Sigma = covariance,
#             group_parameters = group_parameters,
#             joint_parameters = joint_parameters)
# res

## ----eval = FALSE-------------------------------------------------------------
# sres <- summary(res, ci = 0.95)
# sres