## example 1:
# the data consist of ages (in months) at which an
# infant starts to walk alone.
# prepare data
DATA1 <- subset(ZelazoKolb1972, Group != "Control")
# fit unrestricted linear model
fit1.lm <- lm(Age ~ -1 + Group, data = DATA1)
# the variable names can be used to impose constraints on
# the corresponding regression parameters.
coef(fit1.lm)
# constraint syntax: assuming that the walking
# exercises would not have a negative effect of increasing the
# mean age at which a child starts to walk.
myConstraints1 <- ' GroupActive < GroupPassive < GroupNo '
iht(fit1.lm, myConstraints1)
# another way is to first fit the restricted model
fit.restr1 <- restriktor(fit1.lm, constraints = myConstraints1)
iht(fit.restr1)
# Or in matrix notation.
Amat1 <- rbind(c(-1, 0, 1),
c( 0, 1, -1))
myRhs1 <- rep(0L, nrow(Amat1))
myNeq1 <- 0
iht(fit1.lm, constraints = Amat1, rhs = myRhs1, neq = myNeq1)
#########################
## Artificial examples ##
#########################
# generate data
n <- 10
means <- c(1,2,1,3)
nm <- length(means)
group <- as.factor(rep(1:nm, each = n))
y <- rnorm(n * nm, rep(means, each = n))
DATA2 <- data.frame(y, group)
# fit unrestricted linear model
fit2.lm <- lm(y ~ -1 + group, data = DATA2)
coef(fit2.lm)
## example 2: increasing means
myConstraints2 <- ' group1 < group2 < group3 < group4 '
# compute F-test for hypothesis test Type A and compute the tail
# probability based on the parametric bootstrap. We only generate 9
# bootstrap samples in this example; in practice you may wish to
# use a much higher number.
iht(fit2.lm, constraints = myConstraints2, type = "A",
boot = "parametric", R = 9)
# or fit restricted linear model
fit2.con <- restriktor(fit2.lm, constraints = myConstraints2)
iht(fit2.con)
# increasing means in matrix notation.
Amat2 <- rbind(c(-1, 1, 0, 0),
c( 0,-1, 1, 0),
c( 0, 0,-1, 1))
myRhs2 <- rep(0L, nrow(Amat2))
myNeq2 <- 0
iht(fit2.con, constraints = Amat2, rhs = myRhs2, neq = myNeq2,
type = "A", boot = "parametric", R = 9)
## example 3: equality constraints only.
myConstraints3 <- ' group1 = group2 = group3 = group4 '
iht(fit2.lm, constraints = myConstraints3)
# or
fit3.con <- restriktor(fit2.lm, constraints = myConstraints3)
iht(fit3.con)
## example 4:
# combination of equality and inequality constraints.
myConstraints4 <- ' group1 = group2
group3 < group4 '
iht(fit2.lm, constraints = myConstraints4, type = "B", neq.alt = 1)
# fit resticted model and compute model-based bootstrapped
# standard errors. We only generate 9 bootstrap samples in this
# example; in practice you may wish to use a much higher number.
# Note that, a warning message may be thrown because the number of
# bootstrap samples is too low.
fit4.con <- restriktor(fit2.lm, constraints = myConstraints4,
se = "boot.model.based", B = 9)
iht(fit4.con, type = "B", neq.alt = 1)
## example 5:
# restriktor can also be used to define effects using the := operator
# and impose constraints on them. For example, is the
# average effect (AVE) larger than zero?
# generate data
n <- 30
b0 <- 10; b1 = 0.5; b2 = 1; b3 = 1.5
X <- c(rep(c(0), n/2), rep(c(1), n/2))
set.seed(90)
Z <- rnorm(n, 16, 5)
y <- b0 + b1*X + b2*Z + b3*X*Z + rnorm(n, 0, sd = 10)
DATA3 = data.frame(cbind(y, X, Z))
# fit linear model with interaction
fit5.lm <- lm(y ~ X*Z, data = DATA3)
# constraint syntax
myConstraints5 <- ' AVE := X + 16.86137*X.Z;
AVE > 0 '
iht(fit5.lm, constraints = myConstraints5)
# or
fit5.con <- restriktor(fit5.lm, constraints = ' AVE := X + 16.86137*X.Z;
AVE > 0 ')
iht(fit5.con)
# testing equality and/or inequality restrictions in SEM:
#########################
### real data example ###
#########################
# Multiple group path model for facial burns example.
# model syntax with starting values.
burns.model <- 'Selfesteem ~ Age + c(m1, f1)*TBSA + HADS +
start(-.10, -.20)*TBSA
HADS ~ Age + c(m2, f2)*TBSA + RUM +
start(.10, .20)*TBSA '
# constraints syntax
burns.constraints <- 'f2 > 0 ; m1 < 0
m2 > 0 ; f1 < 0
f2 > m2 ; f1 < m1'
# we only generate 2 bootstrap samples in this example; in practice
# you may wish to use a much higher number.
# the double bootstrap was switched off; in practice you probably
# want to set it to "standard".
example6 <- iht(model = burns.model, data = FacialBurns,
R = 2, constraints = burns.constraints,
double.bootstrap = "no", group = "Sex")
example6
##########################
### artificial example ###
##########################
# \donttest{
# Simple ANOVA model with 3 groups (N = 20 per group)
set.seed(1234)
Y <- cbind(c(rnorm(20,0,1), rnorm(20,0.5,1), rnorm(20,1,1)))
grp <- c(rep("1", 20), rep("2", 20), rep("3", 20))
Data <- data.frame(Y, grp)
#create model matrix
fit.lm <- lm(Y ~ grp, data = Data)
mfit <- fit.lm$model
mm <- model.matrix(mfit)
Y <- model.response(mfit)
X <- data.frame(mm[,2:3])
names(X) <- c("d1", "d2")
Data.new <- data.frame(Y, X)
# model
model <- 'Y ~ 1 + a1*d1 + a2*d2'
# fit without constraints
fit <- lavaan::sem(model, data = Data.new)
# constraints syntax: mu1 < mu2 < mu3
constraints <- ' a1 > 0
a1 < a2 '
# we only generate 10 bootstrap samples in this example; in practice
# you may wish to use a much higher number, say > 1000. The double
# bootstrap is not necessary in case of an univariate ANOVA model.
example7 <- iht(model = model, data = Data.new,
start = lavaan::parTable(fit),
R = 10L, double.bootstrap = "no",
constraints = constraints)
example7
# }
Run the code above in your browser using DataLab