library(nlme)
data(Machines)
head(Machines)
## Grouped Data: score ~ Machine | Worker
##   Worker Machine score
## 1      1       A  52.0
## 2      1       A  52.8
## 3      1       A  53.1
## 4      2       A  51.8
## 5      2       A  52.8
## 6      2       A  53.1
machdat = groupedData( score ~ Machine | worker ,
    data = data.frame
        ( 
        score = Machines$score,
        Machine = Machines$Machine,
        worker = Machines$Worker,
        obs=c(rep(1:3,18))),
    labels = list(score = "Score", machine = "Machines"),
    order.groups = FALSE )
mach1 = lme(score ~ Machine-1, data = machdat,
            random = list(worker=pdSymm(~ Machine - 1)), 
            correlation = corSymm(form=~1 | worker/obs), method="ML")
summary(mach1)
## Linear mixed-effects model fit by maximum likelihood
##  Data: machdat 
##        AIC     BIC    logLik
##   241.0082 266.865 -107.5041
## 
## Random effects:
##  Formula: ~Machine - 1 | worker
##  Structure: General positive-definite
##          StdDev    Corr         
## MachineA 3.7174550 MachnA MachnB
## MachineB 7.8707552 0.807        
## MachineC 4.0010768 0.625  0.772 
## Residual 0.9557069              
## 
## Correlation Structure: General
##  Formula: ~1 | worker/obs 
##  Parameter estimate(s):
##  Correlation: 
##   1      2     
## 2 -0.274       
## 3 -0.034  0.045
## Fixed effects: score ~ Machine - 1 
##             Value Std.Error DF  t-value p-value
## MachineA 52.35556  1.578752 46 33.16262       0
## MachineB 60.32222  3.314493 46 18.19953       0
## MachineC 66.27222  1.696696 46 39.05957       0
##  Correlation: 
##          MachnA MachnB
## MachineB 0.794        
## MachineC 0.612  0.763 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.48938841 -0.52730240  0.03136884  0.45026730  2.58923220 
## 
## Number of Observations: 54
## Number of Groups: 6
VarCorr(mach1)
## worker = pdSymm(Machine - 1) 
##          Variance   StdDev    Corr         
## MachineA 13.8194717 3.7174550 MachnA MachnB
## MachineB 61.9487879 7.8707552 0.807        
## MachineC 16.0086154 4.0010768 0.625  0.772 
## Residual  0.9133757 0.9557069
mach2 = lme(score ~ Machine-1, data = machdat,
            random = list(worker=pdSymm(~ Machine - 1)), 
            correlation = corCompSymm(form=~1 | worker/obs), method="ML")
summary(mach2)
## Linear mixed-effects model fit by maximum likelihood
##  Data: machdat 
##        AIC      BIC    logLik
##   237.6593 259.5381 -107.8296
## 
## Random effects:
##  Formula: ~Machine - 1 | worker
##  Structure: General positive-definite
##          StdDev    Corr         
## MachineA 3.7169860 MachnA MachnB
## MachineB 7.8706011 0.806        
## MachineC 4.0005903 0.627  0.774 
## Residual 0.9615755              
## 
## Correlation Structure: Compound symmetry
##  Formula: ~1 | worker/obs 
##  Parameter estimate(s):
##        Rho 
## -0.1367923 
## Fixed effects: score ~ Machine - 1 
##             Value Std.Error DF  t-value p-value
## MachineA 52.35556  1.578767 46 33.16231       0
## MachineB 60.32222  3.314529 46 18.19934       0
## MachineC 66.27222  1.696689 46 39.05974       0
##  Correlation: 
##          MachnA MachnB
## MachineB 0.794        
## MachineC 0.612  0.763 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.45260333 -0.49717339  0.02113936  0.45500017  2.54895024 
## 
## Number of Observations: 54
## Number of Groups: 6
VarCorr(mach2)
## worker = pdSymm(Machine - 1) 
##          Variance   StdDev    Corr         
## MachineA 13.8159851 3.7169860 MachnA MachnB
## MachineB 61.9463622 7.8706011 0.806        
## MachineC 16.0047227 4.0005903 0.627  0.774 
## Residual  0.9246275 0.9615755
mach3 = lme(score ~ Machine-1, data = machdat,
            random = list(worker=pdCompSymm(~ Machine - 1)), 
            correlation = corSymm(form=~1 | worker/obs), method="ML")
summary(mach3)
## Linear mixed-effects model fit by maximum likelihood
##  Data: machdat 
##        AIC      BIC    logLik
##   241.8568 259.7576 -111.9284
## 
## Random effects:
##  Formula: ~Machine - 1 | worker
##  Structure: Compound Symmetry
##          StdDev    Corr       
## MachineA 5.5319899            
## MachineB 5.5319899 0.624      
## MachineC 5.5319899 0.624 0.624
## Residual 0.9557419            
## 
## Correlation Structure: General
##  Formula: ~1 | worker/obs 
##  Parameter estimate(s):
##  Correlation: 
##   1      2     
## 2 -0.275       
## 3 -0.032  0.043
## Fixed effects: score ~ Machine - 1 
##             Value Std.Error DF  t-value p-value
## MachineA 52.35556  2.335433 46 22.41793       0
## MachineB 60.32222  2.335433 46 25.82914       0
## MachineC 66.27222  2.335433 46 28.37685       0
##  Correlation: 
##          MachnA MachnB
## MachineB 0.615        
## MachineC 0.618  0.618 
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.287801061 -0.555270732  0.008721852  0.450551294  2.628294249 
## 
## Number of Observations: 54
## Number of Groups: 6
VarCorr(mach3)
## worker = pdCompSymm(Machine - 1) 
##          Variance   StdDev    Corr       
## MachineA 30.6029124 5.5319899            
## MachineB 30.6029124 5.5319899 0.624      
## MachineC 30.6029124 5.5319899 0.624 0.624
## Residual  0.9134426 0.9557419
mach4 = lme(score ~ Machine-1, data = machdat,
            random = list(worker=pdCompSymm(~ Machine - 1)), 
            correlation = corCompSymm(form=~1 | worker/obs), method="ML")
summary(mach4)
## Linear mixed-effects model fit by maximum likelihood
##  Data: machdat 
##        AIC      BIC    logLik
##   238.5109 252.4338 -112.2555
## 
## Random effects:
##  Formula: ~Machine - 1 | worker
##  Structure: Compound Symmetry
##          StdDev    Corr       
## MachineA 5.5307196            
## MachineB 5.5307196 0.624      
## MachineC 5.5307196 0.624 0.624
## Residual 0.9615753            
## 
## Correlation Structure: Compound symmetry
##  Formula: ~1 | worker/obs 
##  Parameter estimate(s):
##        Rho 
## -0.1367918 
## Fixed effects: score ~ Machine - 1 
##             Value Std.Error DF  t-value p-value
## MachineA 52.35556  2.335043 46 22.42167       0
## MachineB 60.32222  2.335043 46 25.83346       0
## MachineC 66.27222  2.335043 46 28.38159       0
##  Correlation: 
##          MachnA MachnB
## MachineB 0.617        
## MachineC 0.617  0.617 
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.262636022 -0.544930060 -0.006765808  0.446609165  2.573438666 
## 
## Number of Observations: 54
## Number of Groups: 6
VarCorr(mach4)
## worker = pdCompSymm(Machine - 1) 
##          Variance  StdDev    Corr       
## MachineA 30.588860 5.5307196            
## MachineB 30.588860 5.5307196 0.624      
## MachineC 30.588860 5.5307196 0.624 0.624
## Residual  0.924627 0.9615753
mach5 = lme(score ~ Machine-1, data = machdat,
            random = list(worker=pdCompSymm(~ Machine - 1)), 
            correlation = corAR1(form=~1 | worker/obs), method="ML")
summary(mach5)
## Linear mixed-effects model fit by maximum likelihood
##  Data: machdat 
##        AIC      BIC    logLik
##   238.4215 252.3444 -112.2108
## 
## Random effects:
##  Formula: ~Machine - 1 | worker
##  Structure: Compound Symmetry
##          StdDev    Corr       
## MachineA 5.5334801            
## MachineB 5.5334801 0.625      
## MachineC 5.5334801 0.625 0.625
## Residual 0.9614114            
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | worker/obs 
##  Parameter estimate(s):
##        Phi 
## -0.1829829 
## Fixed effects: score ~ Machine - 1 
##             Value Std.Error DF  t-value p-value
## MachineA 52.35556  2.336192 46 22.41063       0
## MachineB 60.32222  2.336192 46 25.82074       0
## MachineC 66.27222  2.336192 46 28.36762       0
##  Correlation: 
##          MachnA MachnB
## MachineB 0.617        
## MachineC 0.619  0.617 
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.264218616 -0.542925969 -0.009085537  0.437344443  2.596943783 
## 
## Number of Observations: 54
## Number of Groups: 6
VarCorr(mach5)
## worker = pdCompSymm(Machine - 1) 
##          Variance   StdDev    Corr       
## MachineA 30.6194017 5.5334801            
## MachineB 30.6194017 5.5334801 0.625      
## MachineC 30.6194017 5.5334801 0.625 0.625
## Residual  0.9243118 0.9614114
# mach6 = lme(score ~ Machine-1, data = machdat,
#             random = list(worker=pdCompSymm(~ Machine - 1)), 
#             weights = varPower(form=~Machine),
#             correlation = corSymm(form=~1 | worker/obs), method="ML")
# summary(mach6)
# VarCorr(mach6)
machA = lme(score ~ Machine-1, data = machdat,
            random = ~ Machine - 1, method="ML")
summary(machA)
## Linear mixed-effects model fit by maximum likelihood
##  Data: machdat 
##        AIC      BIC    logLik
##   236.4178 256.3077 -108.2089
## 
## Random effects:
##  Formula: ~Machine - 1 | worker
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev    Corr         
## MachineA 3.7169532 MachnA MachnB
## MachineB 7.8705146 0.805        
## MachineC 4.0006134 0.625  0.772 
## Residual 0.9615766              
## 
## Fixed effects: score ~ Machine - 1 
##             Value Std.Error DF  t-value p-value
## MachineA 52.35556  1.578753 46 33.16260       0
## MachineB 60.32222  3.314493 46 18.19953       0
## MachineC 66.27222  1.696698 46 39.05952       0
##  Correlation: 
##          MachnA MachnB
## MachineB 0.794        
## MachineC 0.612  0.763 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.40773321 -0.51889887  0.03228994  0.45599259  2.54088421 
## 
## Number of Observations: 54
## Number of Groups: 6
VarCorr(machA)
## worker = pdLogChol(Machine - 1) 
##          Variance   StdDev    Corr         
## MachineA 13.8157407 3.7169532 MachnA MachnB
## MachineB 61.9450000 7.8705146 0.805        
## MachineC 16.0049074 4.0006134 0.625  0.772 
## Residual  0.9246296 0.9615766