# chutepos: commands to estimate variance components for chute data, using 'new' lme
#   multiple lines per parachute

chute2 <- read.table('g:/philip/stat415/chute2.txt',as.is=T,header=T)

names(chute2)		# variable names in the data frame
chute2[1:10,]		# first ten observations

# fit a model with fixed effects: linear regression (intercept and slope)
#                  random effects: intercept and slope, possibly correlated, subject = chute
chute2.lme0 <- lme(strength~position,random=~position|chute,data=chute2)

# fit a model with same fixed effects, but uncorrelated random effects
#   pdDiag specifies a diagonal variance-covariance matrix
chute2.lme1 <- lme(strength~position,random=list(chute=pdDiag(~position)),data=chute2)

# fit a model with only a random slope, -1 suppresses the default intercept
chute2.lme2 <- lme(strength~position,random=~-1+position|chute,data=chute2)

# fit a model with only a random intercept
chute2.lme3 <- lme(strength~position,random=~1|chute,data=chute2)

# see the details for any specific model
summary(chute2.lme0)

# which gives the following output:
summary(chute2.lme0)
Linear mixed-effects model fit by REML
 Data: chute2 
       AIC      BIC    logLik 
  377.1012 386.2594 -182.5506

Random effects:
 Formula:  ~ position | chute
 Structure: General positive-definite
                StdDev   Corr 
(Intercept) 163.063952 (Inter
   position   1.689082 -1    
   Residual  31.887830       

Fixed effects: strength ~ position 
               Value Std.Error DF  t-value p-value 
(Intercept) 1091.744  67.66475 29 16.13461  <.0001
   position    8.962   3.18742 29  2.81165  0.0087
 Correlation: 
         (Intr) 
position -0.37 

Standardized Within-Group Residuals:
       Min         Q1          Med        Q3      Max 
 -1.588375 -0.5861859 -0.002754016 0.3981175 2.177688

Number of Observations: 36
Number of Groups: 6 

# summary model selection and AIC values
anova(chute2.lme3,chute2.lme2,chute2.lme1,chute2.lme0)

# which gives the following output:
 
            Model df      AIC      BIC    logLik   Test  L.Ratio p-value 
chute2.lme3     1  4 373.3452 379.4506 -182.6726                        
chute2.lme2     2  4 424.9654 431.0708 -208.4827                        
chute2.lme1     3  5 375.3452 382.9770 -182.6726 2 vs 3 51.62017  <.0001
chute2.lme0     4  6 377.1012 386.2594 -182.5506 3 vs 4  0.24401  0.6213