Running 'BLUP.tests.R' [170s]
  Running 'OLS.r' [0s]
  Running 'predictionVariance.r' [0s]
  Running 'regress.tests.R' [1s]
  Running 'regressPaper.R' [1s]
Running the tests in 'tests/BLUP.tests.R' failed.
Complete output:
  > set.seed(1001)
  > library(regress)
  > n <- 101
  > x1 <- runif(n)
  > x2 <- seq(0,1,l=n)
  > z1 <- gl(4,10,n)
  > z2 <- gl(6,1,n)
  > 
  > X <- model.matrix(~1 + x1 + x2)
  > Z1 <- model.matrix(~z1-1)
  > Z2 <- model.matrix(~z2-1)
  > 
  > ## Create the individual random and fixed effects
  > beta <- c(1,2,3)
  > eta1 <- rnorm(ncol(Z1),0,10)
  > eta2 <- rnorm(ncol(Z2),0,10)
  > eps <- rnorm(n,0,3)
  > 
  > ## Combine them into a response
  > y <- X %*% beta + Z1 %*% eta1 + Z2 %*% eta2 + eps
  > 
  > ## Fit the same model again
  > model <- regress(y~1 + x1 + x2,~z1 + z2, identity=TRUE,verbose=2)
  1 sigma = 111.2472 111.2472 111.2472 resid llik = -193.5592 
  1 adjusted sigma = 14.95327 14.95327 14.95327
  2 sigma = 48.05741 34.5205 12.87801 resid llik = -185.4081 
  2 adjusted sigma = 43.27607 31.08598 11.59674 delta.llik = 8.151054 
  3 sigma = 70.09679 47.21641 11.03377 resid llik = -184.5278 
  3 adjusted sigma = 69.55275 46.84995 10.94813 delta.llik = 0.8802895 
  4 sigma = 80.45497 53.4053 10.80752 resid llik = -184.4768 
  4 adjusted sigma = 80.4402 53.3955 10.80554 delta.llik = 0.05097438 
  5 sigma = 81.55823 54.06055 10.79318 resid llik = -184.4764 
  5 adjusted sigma = 81.5581 54.06046 10.79316 delta.llik = 0.0004381634 
  > 
  > summary(model)
  Likelihood kernel: K = (Intercept)+x1+x2
  
  Maximized log likelihood with kernel K is  -184.476 
  
  Linear Coefficients:
               Estimate Std. Error
   (Intercept)   -4.573      5.492
   x1             2.477      1.112
   x2             1.967      1.137
  
  Variance Coefficients:
               Estimate Std. Error
            z1   81.671     66.959
            z2   54.127     34.603
            In   10.792      1.609
  
  > names(model)
   [1] "trace"              "llik"               "cycle"             
   [4] "rdf"                "beta"               "beta.cov"          
   [7] "beta.se"            "sigma"              "sigma.cov"         
  [10] "W"                  "Q"                  "fitted"            
  [13] "predicted"          "predictedVariance"  "predictedVariance2"
  [16] "pos"                "Vnames"             "formula"           
  [19] "Vformula"           "Kcolnames"          "model"             
  [22] "Z"                  "X"                  "Sigma"             
  > 
  > b1 <- BLUP(model,RE="z1")
  > 
  > b1$Mean
         z1.1        z1.2        z1.3        z1.4 
    0.3541136   9.0695665   2.9342526 -12.3579327 
  > b1$Variance
      z1.1     z1.2     z1.3     z1.4 
  6.570444 6.570444 6.701733 6.723185 
  > b1$Covariance
           z1.1     z1.2     z1.3     z1.4
  z1.1 6.570444 6.212290 6.200585 6.198651
  z1.2 6.212290 6.570444 6.200585 6.198651
  z1.3 6.200585 6.200585 6.701733 6.185501
  z1.4 6.198651 6.198651 6.185501 6.723185
  > cov2cor(b1$Covariance) ## Large correlation terms
            z1.1      z1.2      z1.3      z1.4
  z1.1 1.0000000 0.9454902 0.9344192 0.9326363
  z1.2 0.9454902 1.0000000 0.9344192 0.9326363
  z1.3 0.9344192 0.9344192 1.0000000 0.9214967
  z1.4 0.9326363 0.9326363 0.9214967 1.0000000
  > 
  > ## Simulate BLUPs given their mean and covariance structure
  > library(MASS)
  > tt <- mvrnorm(n=1000,mu=b1$Mean,Sigma=b1$Covariance)
  > plot(data.frame(tt))
  > matplot(t(tt),type="l",col="grey") ## So there is up and down shift only really.
  > lines(b1$Mean,col="red",lwd=2,type="b")
  > 
  > b2 <- BLUP(model,RE="z2")
  > 
  > b12 <- BLUP(model,RE=c("z1","z2")) ## To examine covariance of z1 and z2 BLUPs
  > 
  > bALL <- BLUP(model) ## To see blups of the measurement error epsilon
  > 
  > ## Plot and compare the random effects and their BLUPs
  > plot(eta1,BLUP(model,RE="z1")$Mean,xlab="Random Effect",ylab="BLUP")
  > plot(eta2,BLUP(model,RE="z2")$Mean,xlab="Random Effect",ylab="BLUP")
  > plot(eps,BLUP(model,RE="In")$Mean,xlab="Random Effect",ylab="BLUP")
  > 
  > 
  > proc.time()
     user  system elapsed 
      0.9     0.1     1.0