Two-stage estimator for non-linear structural equation models
# S3 method for class 'lvmfit'
twostage(
object,
model2,
data = NULL,
predict.fun = NULL,
id1 = NULL,
id2 = NULL,
all = FALSE,
formula = NULL,
std.err = TRUE,
...
)Stage 1 measurement model
Stage 2 SEM
data.frame
Prediction of latent variable
Optional id-variable (stage 1 model)
Optional id-variable (stage 2 model)
If TRUE return additional output (naive estimates)
optional formula specifying non-linear relation
If FALSE calculations of standard errors will be skipped
Additional arguments to lower level functions
m <- lvm(c(x1,x2,x3)~f1,f1~z,
c(y1,y2,y3)~f2,f2~f1+z)
latent(m) <- ~f1+f2
d <- simulate(m,100,p=c("f2,f2"=2,"f1,f1"=0.5),seed=1)
## Full MLE
ee <- estimate(m,d)
## Manual two-stage
if (FALSE) { # \dontrun{
m1 <- lvm(c(x1,x2,x3)~f1,f1~z); latent(m1) <- ~f1
e1 <- estimate(m1,d)
pp1 <- predict(e1,f1~x1+x2+x3)
d$u1 <- pp1[,]
d$u2 <- pp1[,]^2+attr(pp1,"cond.var")[1]
m2 <- lvm(c(y1,y2,y3)~eta,c(y1,eta)~u1+u2+z); latent(m2) <- ~eta
e2 <- estimate(m2,d)
} # }
## Two-stage
m1 <- lvm(c(x1,x2,x3)~f1,f1~z); latent(m1) <- ~f1
m2 <- lvm(c(y1,y2,y3)~eta,c(y1,eta)~u1+u2+z); latent(m2) <- ~eta
pred <- function(mu,var,data,...)
cbind("u1"=mu[,1],"u2"=mu[,1]^2+var[1])
(mm <- twostage(m1,model2=m2,data=d,predict.fun=pred))
#> Estimate Std. Error Z-value P-value
#> Measurements:
#> y2~eta 0.96270 0.12944 7.43758 <1e-12
#> y3~eta 0.97886 0.12738 7.68432 <1e-12
#> Regressions:
#> y1~u1 -0.10923 0.26086 -0.41873 0.6754
#> y1~u2 -0.00916 0.02952 -0.31030 0.7563
#> y1~z -0.09088 0.25976 -0.34985 0.7265
#> eta~u1 1.23462 0.23383 5.27995 1.292e-07
#> eta~u2 0.00912 0.03527 0.25866 0.7959
#> eta~z 0.84531 0.38028 2.22287 0.02622
#> Intercepts:
#> y2 -0.19048 0.18177 -1.04792 0.2947
#> y3 0.00979 0.51253 0.01910 0.9848
#> eta -0.16805 2.38967 -0.07032 0.9439
#> Residual Variances:
#> y1 1.15103 0.24978 4.60814
#> y2 0.97707 0.20206 4.83556
#> y3 1.13661 0.20716 5.48651
#> eta 1.58985 0.41522 3.82896
if (interactive()) {
pf <- function(p) p["eta"]+p["eta~u1"]*u + p["eta~u2"]*u^2
plot(mm,f=pf,data=data.frame(u=seq(-2,2,length.out=100)),lwd=2)
}
## Reduce test timing
## Splines
f <- function(x) cos(2*x)+x+-0.25*x^2
m <- lvm(x1+x2+x3~eta1, y1+y2+y3~eta2, latent=~eta1+eta2)
functional(m, eta2~eta1) <- f
d <- sim(m,500,seed=1,latent=TRUE)
m1 <- lvm(x1+x2+x3~eta1,latent=~eta1)
m2 <- lvm(y1+y2+y3~eta2,latent=~eta2)
mm <- twostage(m1,m2,formula=eta2~eta1,type="spline")
if (interactive()) plot(mm)
nonlinear(m2,type="quadratic") <- eta2~eta1
a <- twostage(m1,m2,data=d)
if (interactive()) plot(a)
kn <- c(-1,0,1)
nonlinear(m2,type="spline",knots=kn) <- eta2~eta1
a <- twostage(m1,m2,data=d)
x <- seq(-3,3,by=0.1)
y <- predict(a, newdata=data.frame(eta1=x))
if (interactive()) {
plot(eta2~eta1, data=d)
lines(x,y, col="red", lwd=5)
p <- estimate(a,f=function(p) predict(a,p=p,newdata=x))$coefmat
plot(eta2~eta1, data=d)
lines(x,p[,1], col="red", lwd=5)
confband(x,lower=p[,3],upper=p[,4],center=p[,1], polygon=TRUE, col=Col(2,0.2))
l1 <- lm(eta2~splines::ns(eta1,knots=kn),data=d)
p1 <- predict(l1,newdata=data.frame(eta1=x),interval="confidence")
lines(x,p1[,1],col="green",lwd=5)
confband(x,lower=p1[,2],upper=p1[,3],center=p1[,1], polygon=TRUE, col=Col(3,0.2))
}
## Reduce test timing
if (FALSE) ## Reduce timing
## Cross-validation example
ma <- lvm(c(x1,x2,x3)~u,latent=~u)
ms <- functional(ma, y~u, value=function(x) -.4*x^2)
#> Error: object 'ma' not found
d <- sim(ms,500)#,seed=1)
#> Error: object 'ms' not found
ea <- estimate(ma,d)
#> Error: object 'ma' not found
mb <- lvm()
mb1 <- nonlinear(mb,type="linear",y~u)
mb2 <- nonlinear(mb,type="quadratic",y~u)
mb3 <- nonlinear(mb,type="spline",knots=c(-3,-1,0,1,3),y~u)
mb4 <- nonlinear(mb,type="spline",knots=c(-3,-2,-1,0,1,2,3),y~u)
ff <- lapply(list(mb1,mb2,mb3,mb4),
function(m) function(data,...) twostage(ma,m,data=data,st.derr=FALSE))
a <- cv(ff,data=d,rep=1)
#> Error in cv(ff, data = d, rep = 1): could not find function "cv"
a
#> Estimate Std. Error Z-value P-value
#> Measurements:
#> y2~eta2 1.03076 0.03983 25.87625 <1e-12
#> y3~eta2 0.99635 0.04147 24.02498 <1e-12
#> Regressions:
#> eta2~eta1_1 2.47344 0.16878 14.65450 <1e-12
#> eta2~eta1_2 -0.48653 0.04469 -10.88630 <1e-12
#> Intercepts:
#> y2 0.02005 0.06576 0.30489 0.7604
#> y3 0.07986 0.06781 1.17768 0.2389
#> eta2 1.16241 0.16835 6.90489 5.024e-12
#> Residual Variances:
#> y1 1.11229 0.10795 10.30394
#> y2 1.00927 0.09809 10.28911
#> y3 1.09169 0.09772 11.17168
#> eta2 1.81941 0.18391 9.89283
# \dontrun{}