library(lava)
m <- lvm()
regression(m) <- c(y1,y2,y3) ~ eta1
regression(m) <- c(z1,z2,z3) ~ eta2
latent(m) <- ~eta1+eta2
regression(m) <- eta2~eta1+x
regression(m) <- eta1~x
labels(m) <- c(eta1=expression(eta[1]),eta2=expression(eta[2]))
plot(m)
d <- sim(m,100,seed=1)
e <- estimate(m,d)
e
## Estimate Std. Error Z-value P-value
## Measurements:
## y2~eta1 0.95462 0.08083 11.80993 <1e-12
## y3~eta1 0.98476 0.08922 11.03722 <1e-12
## z2~eta2 0.97038 0.05368 18.07714 <1e-12
## z3~eta2 0.95608 0.05643 16.94182 <1e-12
## Regressions:
## eta1~x 1.24587 0.11486 10.84694 <1e-12
## eta2~eta1 0.95608 0.18008 5.30910 1.102e-07
## eta2~x 1.11495 0.25228 4.41951 9.893e-06
## Intercepts:
## y2 -0.13896 0.12458 -1.11537 0.2647
## y3 -0.07661 0.13869 -0.55241 0.5807
## eta1 0.15801 0.12780 1.23644 0.2163
## z2 -0.00441 0.14858 -0.02969 0.9763
## z3 -0.15900 0.15731 -1.01076 0.3121
## eta2 -0.14143 0.18380 -0.76949 0.4416
## Residual Variances:
## y1 0.69684 0.14858 4.69004
## y2 0.89804 0.16630 5.40026
## y3 1.22456 0.21182 5.78109
## eta1 0.93620 0.19623 4.77084
## z1 1.41422 0.26259 5.38570
## z2 0.87569 0.19463 4.49934
## z3 1.18155 0.22640 5.21883
## eta2 1.24430 0.28992 4.29195
library(gof)
g <- cumres(e,eta2~eta1)
plot(g)
m <- lvm(c(y1,y2,y3)~u,u~x)
transform(m,u2~u) <- function(x) x^2
regression(m) <- z~u2+u
d <- sim(m,200,p=c("z"=-1,"z~u2"=-0.5),seed=1)
m1 <- lvm(c(y1[0:s],y2[0:s],y3[0:s])~1*u,u~x)
latent(m1) <- ~u
(e1 <- estimate(m1,d))
## Estimate Std. Error Z-value P-value
## Regressions:
## u~x 1.06998 0.08208 13.03542 <1e-12
## Intercepts:
## u -0.08871 0.08753 -1.01344 0.3108
## Residual Variances:
## y1 1.00054 0.07075 14.14214
## u 1.19873 0.15503 7.73233
pp <- function(mu,var,data,...) cbind(u=mu[,"u"],u2=mu[,"u"]^2+var["u","u"])
(e <- measurement.error(e1, z~1+x, data=d, predictfun=pp))
## Estimate Std.Err 2.5% 97.5% P-value
## (Intercept) -1.1181 0.13795 -1.3885 -0.8477 5.273e-16
## x -0.0537 0.13213 -0.3127 0.2053 6.844e-01
## u 1.0039 0.11504 0.7785 1.2294 2.609e-18
## u2 -0.4718 0.05213 -0.5740 -0.3697 1.410e-19
f <- function(p) p[1]+p["u"]*u+p["u2"]*u^2
u <- seq(-1,1,length.out=100)
plot(e, f, data=data.frame(u))
m <- lvm(y~x,c~1)
regression(m) <- c(y,x)~z
eventTime(m) <- t~min(y=1,c=0)
transform(m,S~t+status) <- function(x) survival::Surv(x[,1],x[,2])
plot(m)
onerun <- function(...) {
d <- sim(m,100)
m0 <- lvm(S~x+z,x~z)
e <- estimate(m0,d,estimator="glm")
vec(summary(effects(e,S~z))$coef[,1:2])
}
val <- sim(onerun,100)
summary(val, estimate=1:4, se=5:8, short=TRUE)
## 100 replications Time: 8.98s
##
## Total.Estimate Direct.Estimate Indirect.Estimate S~x~z.Estimate
## Mean 2.00547 1.00826 0.99721 0.99721
## SD 0.17529 0.17928 0.18200 0.18200
## SE 0.18540 0.17935 0.16539 0.16539
## SE/SD 1.05767 1.00041 0.90876 0.90876
##
## Min 1.63703 0.57336 0.60752 0.60752
## 2.5% 1.70218 0.64867 0.66202 0.66202
## 50% 1.97280 1.03902 0.98842 0.98842
## 97.5% 2.38755 1.39713 1.38015 1.38015
## Max 2.42266 1.45465 1.52908 1.52908
##
## Missing 0.00000 0.00000 0.00000 0.00000
val <- sim(val,500) ## Add 500 simulations
plot(val,estimate=c("Total.Estimate","Indirect.Estimate"),
true=c(2,1),se=c("Total.Std.Err","Indirect.Std.Err"))
