R Under development (unstable) (2025-05-08 r88190 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> ## Test 1, Random Effects Model
> 
> library(nlme)
> library(regress)
> data(Oats)
> names(Oats) <- c("B","V","N","Y")
> Oats$N <- as.factor(Oats$N)
> 
> ## Using regress
> oats.reg <- regress(Y~N+V,~B+I(B:V),identity=TRUE,verbose=1,data=Oats)
1 sigma = 732.1964 732.1964 732.1964 resid llik = -215.1927 
1 adjusted sigma = 157.9849 157.9849 157.9849
2 sigma = 186.231 133.8389 160.2718 resid llik = -215.0362 
2 adjusted sigma = 185.6373 133.4123 159.761 delta.llik = 0.1565299 
3 sigma = 205.8252 116.8087 161.7195 resid llik = -214.981 
3 adjusted sigma = 205.6644 116.7175 161.5931 delta.llik = 0.05518008 
4 sigma = 213.5958 110.3954 162.4623 resid llik = -214.975 
4 adjusted sigma = 213.5887 110.3917 162.4569 delta.llik = 0.005981576 
5 sigma = 214.3882 109.7628 162.5486 resid llik = -214.9749 
5 adjusted sigma = 214.3882 109.7628 162.5486 delta.llik = 6.213704e-05 
> summary(oats.reg)
Likelihood kernel: K = (Intercept)+N0.2+N0.4+N0.6+VMarvellous+VVictory

Maximized log likelihood with kernel K is  -214.975 

Linear Coefficients:
             Estimate Std. Error
 (Intercept)   79.917      8.220
 N0.2          19.500      4.250
 N0.4          34.833      4.250
 N0.6          44.000      4.250
 VMarvellous    5.292      7.079
 VVictory      -6.875      7.079

Variance Coefficients:
             Estimate Std. Error
      B       214.468    168.794
      I(B:V)  109.700     67.741
      In      162.558     32.189

> 
> ## Using lme
> oats.lme <- lme(Y~N+V,random=~1|B/V,data=Oats,method="REML")
> summary(oats.lme)
Linear mixed-effects model fit by REML
  Data: Oats 
       AIC      BIC    logLik
  586.0688 605.7756 -284.0344

Random effects:
 Formula: ~1 | B
        (Intercept)
StdDev:    14.64487

 Formula: ~1 | V %in% B
        (Intercept) Residual
StdDev:    10.47345 12.74986

Fixed effects:  Y ~ N + V 
               Value Std.Error DF   t-value p-value
(Intercept) 79.91667  8.220350 51  9.721808  0.0000
N0.2        19.50000  4.249955 51  4.588284  0.0000
N0.4        34.83333  4.249955 51  8.196166  0.0000
N0.6        44.00000  4.249955 51 10.353051  0.0000
VMarvellous  5.29167  7.078912 10  0.747525  0.4720
VVictory    -6.87500  7.078912 10 -0.971194  0.3544
 Correlation: 
            (Intr) N0.2   N0.4   N0.6   VMrvll
N0.2        -0.259                            
N0.4        -0.259  0.500                     
N0.6        -0.259  0.500  0.500              
VMarvellous -0.431  0.000  0.000  0.000       
VVictory    -0.431  0.000  0.000  0.000  0.500

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8413413 -0.6627975 -0.0669429  0.6382246  1.6606681 

Number of Observations: 72
Number of Groups: 
       B V %in% B 
       6       18 
> 
> if((oats.lme$sigma^2 - oats.reg$sigma[3])^2>0.0001) stop("Error 1 - doesn't match lme")
> 
> b1 <- ranef(oats.lme)
> b2 <- BLUP(oats.reg)
> if(sum((unlist(b1) - b2$Mean)^2)>0.0002) stop("Error with BLUP")
> 
> ## Test 2, Multivariate Model
> 
> library(regress)
> library(MASS)
> n <- 100
> mu <- c(1,2)
> sigma <- c(10,10,5)
> Sigma <- matrix(c(sigma[1],sigma[3],sigma[3],sigma[2]),2,2)
> Y <- mvrnorm(n,mu,Sigma)
> 
> y <- c(Y[,1],Y[,2])
> X <- kronecker(diag(1,2),rep(1,n))
> 
> V1 <- matrix(c(1,0,0,0),2,2)
> V2 <- matrix(c(0,0,0,1),2,2)
> V3 <- matrix(c(0,1,1,0),2,2)
> 
> sig1 <- kronecker(V1,diag(1,n))
> sig2 <- kronecker(V2,diag(1,n))
> gam <- kronecker(V3,diag(1,n))
> 
> reg.obj <- regress(y~X-1,~sig1+sig2+gam,identity=FALSE,verbose=1,start=c(10,10,5))
1 sigma = 10 10 5 resid llik = -302.2268 
1 adjusted sigma = 8.994715 8.994715 4.497357
2 sigma = 9.001019 7.978671 3.487618 resid llik = -301.1933 
2 adjusted sigma = 8.983818 7.963423 3.480953 delta.llik = 1.033495 
3 sigma = 8.99993 7.875542 3.385977 resid llik = -301.1827 
3 adjusted sigma = 8.999754 7.875388 3.385911 delta.llik = 0.01062619 
4 sigma = 9.001523 7.866738 3.376473 resid llik = -301.1826 
4 adjusted sigma = 9.001522 7.866737 3.376473 delta.llik = 0.0001062922 
> summary(reg.obj)
Likelihood kernel: K = X1+X2

Maximized log likelihood with kernel K is  -301.183 

Linear Coefficients:
      Estimate Std. Error
   X1    0.543       0.30
   X2    1.761       0.28

Variance Coefficients:
      Estimate Std. Error
 sig1    9.002      1.279
 sig2    7.866      1.118
 gam     3.376      0.911

> 
> proc.time()
   user  system elapsed 
   0.53    0.10    0.62