unit tests updated

.gitignore

@@ -34,3 +34,4 @@ test
 *.dvi
 _region*
 symbols.rds
+Rplots.pdf

.travis.yml

@@ -27,6 +27,11 @@ r_binary_packages:
  - geepack
  - testthat
  - png
+ - MatrixModels
+ - SparseM
+ - lazyeval
+ - quantreg
+ - zoo
 
 r_github_packages:
  - jimhester/covr

tests/testthat/test-constrain.R

@@ -36,18 +36,14 @@ test_that("constrain (Fishers z-transform)",{
 })
 
 
-test_that("Multiple group constraints I", {
-    set.seed(1)
-    m1 <- lvm(y[m:v] ~ f(x,beta)+f(z,beta2))
-    d1 <- sim(m1,500); d2 <- sim(m1,500)
-    constrain(m1,beta2~psi) <- function(x) 2*x
-    m2 <- lvm(y[m:v] ~ f(x,beta2) + z)
-    constrain(m2,beta2~psi) <- function(x) 2*x
-    mg <- multigroup(list(m1,m2),list(d1,d2))
-    ee <- estimate(mg)
-    expect_true(length(coef(ee))==5)    
-    expect_equivalent(constraints(ee)[1],2*coef(ee)["1@psi"]) # Est
-    expect_equivalent(constraints(ee)[2],2*coef(ee,2)[[1]]["psi",2]) # Std.Err    
+test_that("Non-linear in exogenous variables", {
+    d <- data.frame(x=1:5,y=c(0.5,1.5,2,3,3.5))
+    m <- lvm(y[m] ~ 1)
+    addvar(m) <- ~x
+    parameter(m) <- ~a+b
+    constrain(m,m~a+b+x) <- function(z) z[1]+z[2]*z[3]
+    e <- estimate(m,d,control=list(method="NR"))
+    expect_true(mean(coef(lm(y~x,d))-coef(e)[c("a","b")])^2<1e-4)
 })
 
 
@@ -66,13 +62,28 @@ test_that("Probit constraints", {
     }   
 })
 
+
+test_that("Multiple group constraints I", {
+    set.seed(1)
+    m1 <- lvm(y[m:v] ~ f(x,beta)+f(z,beta2))
+    d1 <- sim(m1,500); d2 <- sim(m1,500)
+    constrain(m1,beta2~psi) <- function(x) 2*x
+    m2 <- lvm(y[m:v] ~ f(x,beta2) + z)
+    constrain(m2,beta2~psi) <- function(x) 2*x
+    mg <- multigroup(list(m1,m2),list(d1,d2))
+    ee <- estimate(mg)
+    expect_true(length(coef(ee))==5)    
+    expect_equivalent(constraints(ee)[1],2*coef(ee)["1@psi"]) # Est
+    expect_equivalent(constraints(ee)[2],2*coef(ee,2)[[1]]["psi",2]) # Std.Err    
+})
+
 test_that("Multiple group constraints II", {
   data(twindata)
   twinwide <- reshape(twindata,direction="wide",
                       idvar="id",timevar="twinnum")
   l <- lvm(~bw.1+bw.2)
   covariance(l) <- bw.1 ~ bw.2
-  e <- estimate(l,subset(twinwide,zyg.1=="MZ"))
+  e <- estimate(l,subset(twinwide,zyg.1=="MZ"),control=list(method="NR"))
   B <- cbind(1,-1); colnames(B) <- c("bw.1,bw.1","bw.2,bw.2")
   lava::compare(e,contrast=B)
   B2 <- rbind(c(1,-1,0,0),c(0,0,1,-1))
@@ -89,22 +100,19 @@ test_that("Multiple group constraints II", {
 
   DZ <- subset(twinwide,zyg.1=="MZ")
   MZ <- subset(twinwide,zyg.1=="DZ")
-  e <- estimate(l,MZ)
-  e2 <- estimate(l2,DZ)
+  ## e <- estimate(l,MZ)
+  ## e2 <- estimate(l2,DZ)
 
   parameter(l) <- ~r2
   parameter(l2) <- ~r1
-  ee <- estimate(list(MZ=l,DZ=l2),list(MZ,DZ))
-  constrain(ee,h~r1+r2) <- function(x) 2*(x[1]-x[2])
+  ee <- estimate(list(MZ=l,DZ=l2),list(MZ,DZ),control=list(method="NR"))
+  constrain(ee,h~r2+r1) <- function(x) 2*(x[1]-x[2])
   ce <- constraints(ee)
   expect_true(length(coef(ee))==4)
   expect_true(nrow(ce)==1)
   expect_true(all(!is.na(ce)))
   expect_true(mean(score(ee)^2)<1e-4)
 })
 
-## test_that("text",{
-##   ##  expect_output(g,"p=12")
-## })
 
 

tests/testthat/test-inference.R

@@ -88,8 +88,8 @@ test_that("equivalence", {
     regression(m) <- y2~x
     e <- estimate(m,d)
     ##eq <- equivalence(e,y1~x,k=1)
-    eq <- equivalence(e,y2~x,k=1)
-    print(eq)
+    suppressMessages(eq <- equivalence(e,y2~x,k=1))
+    expect_output(print(eq),"y1,y2")
     expect_true(all(c("y1","y3")%in%eq$equiv[[1]][1,]))
 })
 
@@ -101,8 +101,6 @@ test_that("multiple testing", {
     expect_equivalent(ci1[,1],ci2[,1])
     expect_true(all(ci1[,2]<ci2[,2]))
     expect_true(all(ci1[,3]>ci2[,3]))
-
-
 })
 
 
@@ -131,8 +129,8 @@ test_that("Bootstrap", {
     y <- rep(c(0,1),each=5)
     x <- 1:10
     e <- estimate(y~x)
-    B1 <- bootstrap(e,R=2,silent=TRUE)
-    B2 <- bootstrap(e,R=2,silent=TRUE,bollenstine=TRUE)
+    B1 <- bootstrap(e,R=2,silent=TRUE,mc.cores=1)
+    B2 <- bootstrap(e,R=2,silent=TRUE,bollenstine=TRUE,mc.cores=1)
     expect_false(B1$bollenstine)
     expect_true(B2$bollenstine)
     expect_true(nrow(B1$coef)==2)
@@ -202,7 +200,7 @@ test_that("zero-inflated binomial regression (zib)", {
 
 test_that("Optimization", {
     m <- lvm(y~x+z)
-    d <- simulate(m,10,seed=1)
+    d <- simulate(m,20,seed=1)
     e1 <- estimate(m,d,control=list(method="nlminb0"))
     e2 <- estimate(m,d,control=list(method="NR"))
     expect_equivalent(round(coef(e1),3),round(coef(e2),3))
@@ -260,3 +258,49 @@ test_that("Prediction, random intercept", {
     expect_true(mse(u01,u02)<1e-9)    
 })
 
+
+
+test_that("Random slope model", {
+
+    set.seed(1)
+    m <- lvm()
+    regression(m) <- y1 ~ 1*u+1*s
+    regression(m) <- y2 ~ 1*u+2*s
+    regression(m) <- y3 ~ 1*u+3*s
+    latent(m) <- ~u+s
+    d <- sim(m,20)
+
+    dd <- mets::fast.reshape(d)
+    library(lme4)
+    l <- lmer(y~ 1+num +(1+num|id),dd,REML=FALSE)
+    sl <- lava:::varcomp(l,profile=FALSE)
+
+    d0 <- mets::fast.reshape(dd,id="id")
+    m0 <- lvm(c(y1[0:v],y2[0:v],y3[0:v])~1*u)
+    addvar(m0) <- ~num1+num2+num3
+    covariance(m0) <- u~s
+    latent(m0) <- ~s+u
+    regression(m0) <- y1 ~ num1*s
+    regression(m0) <- y2 ~ num2*s
+    regression(m0) <- y3 ~ num3*s
+    system.time(e0 <- estimate(m0,d0,param="none",control=list(trace=0)))
+    system.time(e1 <- estimate(m1,d0,param="none"))
+
+    expect_true(mean(sl$coef-coef(e0)[c("u","s")])^2<1e-5)
+    expect_true((logLik(l)-logLik(e0))^2<1e-5)    
+    expect_true(mean(diag(sl$varcomp)-coef(e0)[c("u,u","s,s")])^2<1e-5)
+
+    m1 <- lvm(c(y1[0:v],y2[0:v],y3[0:v])~1*u)
+    addvar(m1) <- ~num1+num2+num3
+    covariance(m1) <- u~s
+    latent(m1) <- ~s+u
+    regression(m1) <- y1 ~ b1*s
+    regression(m1) <- y2 ~ b2*s
+    regression(m1) <- y3 ~ b3*s
+    constrain(m1,b1~num1) <- function(x) x
+    constrain(m1,b2~num2) <- function(x) x
+    constrain(m1,b3~num3) <- function(x) x
+    system.time(e1 <- estimate(m1,d0,param="none",p=coef(e0)))
+    expect_true((logLik(e1)-logLik(e0))^2<1e-5)    
+
+})

tests/testthat/test-multigroup.R

@@ -54,6 +54,41 @@ test_that("Multiple group, missing data analysis", {
 })
 
 
+test_that("Multiple group, constraints", {
+    m1 <- lvm(y ~ f(x,beta)+f(z,beta2))
+    m2 <- lvm(y ~ f(x,psi) + z)
+    ## And simulate data from them
+    set.seed(1)
+    d1 <- sim(m1,100)
+    d2 <- sim(m2,100)
+    ## Add 'non'-linear parameter constraint
+    constrain(m2,psi ~ beta2) <- function(x) x
+    ## Add parameter beta2 to model 2, now beta2 exists in both models
+    parameter(m2) <- ~ beta2
+    ee <- estimate(list(m1,m2),list(d1,d2))
+
+    m <- lvm(y1 ~ x1 + beta2*z1)
+    regression(m) <- y2 ~ beta2*x2 + z2
+    d <- cbind(d1,d2); names(d) <- c(paste0(names(d1),1),paste0(names(d1),2))
+    e <- estimate(m,d)
+
+    b1 <- coef(e,2)["beta2",1]
+    b2 <- constraints(ee)[1]
+    expect_true(mean((b1-b2)^2)<1e-5)
+
+    ## "Multiple group, constraints (non-linear in x)
+    m <- lvm(y[m:v] ~ 1)
+    addvar(m) <- ~x
+    parameter(m) <- ~a+b
+    constrain(m,m~a+b+x) <- function(z) z[1]+z[2]*z[3]
+    ee <- estimate(list(m,m),list(d1[1:5,],d1[6:10,]))
+    b1 <- coef(lm(y~x,d1[1:10,]))
+    b2 <- coef(ee)[c("1@a","1@b")]
+    expect_true(mean(b1-b2)^2<1e-4)
+
+})
+
+
 
 
 

tests/testthat/test-plot.R

@@ -37,21 +37,6 @@ if (requireNamespace("visualTest") && requireNamespace("png")) {
     }
 
 
-    d1 <- gropen()    
-    plot(1,xlim=c(-2,2),ylim=c(-2,2))
-    points(0.5,0.5)
-    points(1,1)
-    dev.off()
-
-    d2 <- gropen()
-    plot(1,xlim=c(-2,2),ylim=c(-2,2))
-    points(1,1)
-    points(0.5,0.5)
-    dev.off()
-
-    grcompare(d1,d2)
-
-
     test_that("plotConf", {
         set.seed(1)
         x <- rnorm(50)
@@ -99,11 +84,34 @@ if (requireNamespace("visualTest") && requireNamespace("png")) {
         dev.off()
 
         expect_true(grcompare(d1,d2,threshold=10))
-
-        ##forestplot(coef(l0),confint(l0)[,1],confint(l0)[,2],axes=FALSE)
-
+
   })
 
 
+    test_that("forestplot", {
+        set.seed(1)
+        K <- 20
+        est <- rnorm(K); est[c(3:4,10:12)] <- NA        
+        se <- runif(K,0.2,0.4)        
+        x <- cbind(est,est-2*se,est+2*se,runif(K,0.5,2))        
+        rownames(x) <- unlist(lapply(letters[seq(K)],function(x) paste(rep(x,4),collapse="")))
+        rownames(x)[which(is.na(est))] <- ""
+        signif <- sign(x[,2])==sign(x[,3])        
+        forestplot(x)
+        ## TODO
+    })
+
+    test_that("plot.sim", {
+        onerun2 <- function(a,b,...) {
+            return(cbind(a=a,b=b,c=a-1,d=a+1))
+        }
+        R <- data.frame(a=1:2,b=3:4)
+        val2 <- sim(onerun2,R=R,type=0)
+        plot(val2)
+        plot(val2,plot.type="single")
+        density(val2)
+        ## TODO
+    })
+
 
 }

tests/testthat/test-sim.R

@@ -12,7 +12,7 @@ test_that("Missing", {
 })
 
 
-test_that("sim.default", {
+test_that("sim.default I", {
     m <- lvm(y~x+e)
     distribution(m,~y) <- 0
     distribution(m,~x) <- uniform.lvm(a=-1.1,b=1.1)
@@ -31,13 +31,33 @@ test_that("sim.default", {
 
     val <- sim(onerun,R=2,b0=1,n=10,messages=0,mc.cores=1)
     expect_true(nrow(val)==2)
-    val <- sim(val,R=2,b0=1,n=10,mc.cores=1,messages=0) ## append results
+    val <- sim(val,R=2,b0=1,n=10,type=0) ## append results
     expect_true(nrow(val)==4)
 
     s1 <- summary(val,estimate=c(1,1),confint=c(3,4,6,7),true=c(1,1),names=c("Model","Sandwich"))
     expect_true(length(grep("Coverage",rownames(s1)))>0)
     expect_equivalent(colnames(s1),c("Model","Sandwich"))
 
+    val <- sim(onerun,R=2,cl=TRUE,seed=1,messages=0)
+    expect_true(val[1,1]!=val[1,2])
+
+    onerun2 <- function(a,b,...) {
+        return(cbind(a=a,b=b,c=a-1,d=a+1))
+    }
+    R <- data.frame(a=1:2,b=3:4)
+    val2 <- sim(onerun2,R=R,messages=1,mc.cores=2)
+    expect_true(all(R-val2[,1:2]==0))
+    res <- summary(val2)
+    expect_equivalent(res["Mean",],c(1.5,3.5,0.5,2.5))
+
+    expect_output(print(val2[1,]),"a b c d")
+    expect_output(print(val2[1,]),"1 3 0 2")
+
+    res <- summary(val2,estimate="a",se="b",true=1.5,confint=c("c","d"))
+    expect_true(res["Coverage",]==1)
+    expect_true(res["SE/SD",]==mean(val2[,"b"])/sd(val2[,"a"]))
+
 })
 
 
+