Nothing
## ----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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.