Carry-over balanced designs for 8 treatments

Those are Williams designs you might say, but it has become clear to me that Williams designs are just a subset of all carry-over balanced designs. Not through hard work of mine, comments by Apn on my previous post Creating Williams designs with even number of products lead to this. Among other things, Apn claimed there were designs which did not confirm to the well known symmetry that each row is the reflection of another row. Related to that was a claim that there is a nine treatment square carry-over balanced design. I made my algorithm faster and in this post I can confirm the former, but the latter is too long a calculation for the algorithm I used.

Designs

There are 12 designs, as given below. Numbers 4 and 5 can be obtained from each other by certain recoding of the treatments. The last four are the non-symmetrical ones.
[[1]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    4    1    6    3    8    5    7
[3,]    3    1    5    2    7    4    8    6
[4,]    4    6    2    8    1    7    3    5
[5,]    5    3    7    1    8    2    6    4
[6,]    6    8    4    7    2    5    1    3
[7,]    7    5    8    3    6    1    4    2
[8,]    8    7    6    5    4    3    2    1

[[2]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    4    8    3    6    1    5    7
[3,]    3    1    4    7    2    5    8    6
[4,]    4    6    2    8    1    7    3    5
[5,]    5    3    7    1    8    2    6    4
[6,]    6    8    5    2    7    4    1    3
[7,]    7    5    1    6    3    8    4    2
[8,]    8    7    6    5    4    3    2    1

[[3]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    4    1    3    6    8    5    7
[3,]    3    8    4    7    2    5    1    6
[4,]    4    6    2    8    1    7    3    5
[5,]    5    3    7    1    8    2    6    4
[6,]    6    1    5    2    7    4    8    3
[7,]    7    5    8    6    3    1    4    2
[8,]    8    7    6    5    4    3    2    1

[[4]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    4    8    6    3    1    5    7
[3,]    3    8    5    2    7    4    1    6
[4,]    4    6    2    8    1    7    3    5
[5,]    5    3    7    1    8    2    6    4
[6,]    6    1    4    7    2    5    8    3
[7,]    7    5    1    3    6    8    4    2
[8,]    8    7    6    5    4    3    2    1

[[5]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    5    1    6    3    8    4    7
[3,]    3    1    4    2    7    5    8    6
[4,]    4    6    2    8    1    7    3    5
[5,]    5    3    7    1    8    2    6    4
[6,]    6    8    5    7    2    4    1    3
[7,]    7    4    8    3    6    1    5    2
[8,]    8    7    6    5    4    3    2    1

[[6]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    5    8    3    6    1    4    7
[3,]    3    1    5    7    2    4    8    6
[4,]    4    6    2    8    1    7    3    5
[5,]    5    3    7    1    8    2    6    4
[6,]    6    8    4    2    7    5    1    3
[7,]    7    4    1    6    3    8    5    2
[8,]    8    7    6    5    4    3    2    1

[[7]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    5    1    3    6    8    4    7
[3,]    3    8    5    7    2    4    1    6
[4,]    4    6    2    8    1    7    3    5
[5,]    5    3    7    1    8    2    6    4
[6,]    6    1    4    2    7    5    8    3
[7,]    7    4    8    6    3    1    5    2
[8,]    8    7    6    5    4    3    2    1

[[8]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    5    8    6    3    1    4    7
[3,]    3    8    4    2    7    5    1    6
[4,]    4    6    2    8    1    7    3    5
[5,]    5    3    7    1    8    2    6    4
[6,]    6    1    5    7    2    4    8    3
[7,]    7    4    1    3    6    8    5    2
[8,]    8    7    6    5    4    3    2    1

[[9]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    7    1    8    3    5    4    6
[3,]    3    1    5    7    6    8    2    4
[4,]    4    8    7    5    2    1    6    3
[5,]    5    3    6    2    8    4    1    7
[6,]    6    5    8    1    4    7    3    2
[7,]    7    4    2    6    1    3    8    5
[8,]    8    6    4    3    7    2    5    1

[[10]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    7    6    3    8    1    5    4
[3,]    3    6    8    5    2    4    1    7
[4,]    4    3    5    7    1    8    6    2
[5,]    5    8    2    1    3    7    4    6
[6,]    6    1    4    8    7    3    2    5
[7,]    7    5    1    6    4    2    8    3
[8,]    8    4    7    2    6    5    3    1

[[11]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    8    5    1    7    3    6    4
[3,]    3    5    4    6    1    8    2    7
[4,]    4    1    6    8    3    7    5    2
[5,]    5    7    1    3    2    4    8    6
[6,]    6    3    8    7    4    2    1    5
[7,]    7    6    2    5    8    1    4    3
[8,]    8    4    7    2    6    5    3    1

[[12]]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    2    3    4    5    6    7    8
[2,]    2    8    5    7    4    1    3    6
[3,]    3    5    2    6    8    7    1    4
[4,]    4    7    6    2    1    5    8    3
[5,]    5    4    8    1    6    3    2    7
[6,]    6    1    7    5    3    8    4    2
[7,]    7    3    1    8    2    4    6    5
[8,]    8    6    4    3    7    2    5    1

Calculation speed

Up to 5 treatments 100 repeats are done in the blink of an eye. For 6 it is doable (0.025 sec per run). For 7 treatments I did 10 repeats (1.6 sec per run), while for 8 one run took almost 800 sec. Clearly 9 treatments is too much.
> microbenchmark(gendesign(3),gendesign(4),gendesign(5),gendesign(6))
Note: no visible global function definition for 'error'
Unit: microseconds
         expr       min         lq    median         uq       max neval
 gendesign(3)   107.029   113.2605   119.492   125.3565   240.449   100
 gendesign(4)   343.081   355.5430   369.472   405.7585   544.677   100
 gendesign(5)  1141.403  1163.0290  1202.615  1258.6955  1472.755   100
 gendesign(6) 24753.851 25352.7755 25595.058 26917.8980 35721.433   100
> system.time(gendesign(7))
   user  system elapsed
   1.64    0.00    1.67
> microbenchmark(gendesign(7),times=10L)
Unit: seconds
         expr      min       lq   median       uq      max neval
 gendesign(7) 1.624484 1.647051 1.657468 1.669608 1.768613    10
> system.time(gendesign(8))
   user  system elapsed
 793.72    0.17  805.82

Code

Code is below. Running the JIT saves half of the time. It is basically a simplified version of my old algorithm, but with a lot of small modifications, basically avoiding loops, ifs and a load of trial and error for speed.

library(microbenchmark)

nextpos <- function(desmat) which(desmat==0,arr.ind=TRUE)

gendesign <- function(n=6) {

  nr <- as.integer(n)

  nc <- nr

  

  desmat <- matrix(0L,nrow=nr,ncol=nc)

  desmat[1,] <- 1L:nc

  desmat[,1] <- 1L:nr

  

  carover <- matrix(TRUE,nrow=nr,ncol=nc)

  for (i in 1L:(nc-1L))  carover[i,i+1] <- FALSE

  todo <- nextpos(desmat)

  

  desobject <- list(desmat=desmat,carover = carover,nc=1L:n,n=n,

      index =1L,npos=nrow(todo),

      row=todo[,1L],

      col=todo[,2L])

  desresult <- list()

  addpoint(desobject,desresult)

}

modify <- function(desobject,row,col,i,previous) {

  desobject$desmat[row,col] <- i

  desobject$carover[previous,i] <- FALSE

  desobject$index <- desobject$index + 1L

  desobject}

addpoint <- function(desobject,desresult) {

  if (desobject$index>desobject$npos) {

    l <- length(desresult)

    desresult[[l+1]] <- desobject$desmat

  #  cat('#')

    return(desresult)

  } 

  row <- desobject$row[desobject$index]

  col <- desobject$col[desobject$index]

  previous <- desobject$desmat[row,col-1L]

  avoid <- c(desobject$desmat[row,],

      desobject$desmat[,col])

  nc <- desobject$nc[!is.element(desobject$nc,avoid) ]

  nc <- nc[desobject$carover[previous,nc]]

  for (i in nc) {

       desresult <- addpoint(modify(desobject,row,col,i,previous)

       ,desresult)

  }

  desresult

}

library(compiler)

enableJIT(3)

Prices of houses in the Netherlands

The last couple of days I read a number of times about stabilization in house prices which had been dropping due to the crisis. And you get hit by numbers such as change against Q2 2013 or Q3 2012. These are accompanied by reasons why this or that quarter may be special so changes may be off. To be honest, after years of crisis I for one have no clue how prices have changed overall. Hence this small good news was a good reason to grab some data and look for myself.

Data

The data is obtained from the website of NVM, the largest Dutch organization of real estate agents. Their website contains loads of reports and data on the housing market, most of them in Dutch. Since I wanted analyze data I started with grabbing historical data. These are house prices since 1985. Note that the good news was that the number of transactions seemed to increase, not the prices.

The regions are chosen because they reflect some specific properties. These descriptions are purely my own.
Amsterdam is the financial heart and may be sensitive to the financial markets.
Den Haag (The Hague) is the government and related institutions
Rotterdam the harbor, with logistics, may be sensitive to neighbor Germany
Waterland typical location if you work in Amsterdam but don't want to live there
Almere another typical location to work in Amsterdam, a new grow city, almost doubled in size since 1995
Zuidwest Drente is the countryside
Ede eo is food valley, the more high tech food companies may be found there
Zeeuwse eilanden is more tourism and fishing
Eindhoven is a region with more high tech industries

The historical data is a .pdf per NVM region, covers 1985-2012. These are a bit strange files. I was unable to copy paste the second page of the files, at least into LibreOffice Calc. Plan B was to convert via pdftotext which can be found as part of xpdf. The resulting .txt file is no beauty either, with loads of data on one line, separated by some keywords, then data every other line.
Data for 2013 is not historical. Q3 is in a small report for each region, the new data, which luckily also contains Q2. Q1 was extracted by getting same report for Q2. It seemed most easy to manually extract these in a table and then read them in R.
The code for reading the data is in the appendix.

Analysis

Raw data 

The first plot contains all data, on a logarithmic scale for price. From this plot it seems the prices have increased gigantic. From just over €50000 to €250000 in 30 years. Its also seems there is a linear phase (hence exponential growth) from say 1985 to 1995 followed by an even faster increase till 2002, which flattened till 2007 after which the decrease set in. A long term interpretation might be that the last couple of years prices corrected after the crazy times around the turn of the century.

library(ggplot2)

p <- ggplot(all, aes(x=time, y=price))

p + geom_line() +

    scale_y_log10(limits=c(50,300),breaks=seq(50,300,by=50)) + 

    scale_x_continuous(breaks=seq(1990,2010,by=10),'Year') +

    facet_wrap(~region)

Smoothed data

I love my smoothers, so using them is almost second nature. I also took the opportunity to restrict myself to the year 2000 to 2013 and moved to a linear scale. These plots do suggest that Amsterdam, Den Haag and Eindhoven found the way up again, but there are just as many regions which keep on dropping in price.

p + stat_smooth(method='loess',span=.2) +

    scale_y_continuous(limits=c(130,300),

        'Price (thousands Euros)') + 

    scale_x_continuous(limits=c(2000,2013.8),'Year') +

    facet_wrap(~region)

Forecasting

I have seen many a blog with forecast, all the more reason to try to use it. It was easy to use. However, I have some fears for anything 'auto'. It hides a lot. Further analysis (not shown) showed many competing models and models which could not be estimated. Some of these models are close in aicc, yet predict differently. However, time series has been a bit too long ago to really believe I can do better within a reasonable time. Good thing about 'auto' is my personal desire to obtain certain effects is also out of the picture.

The plots show an increases in Den Haag, maybe in Amsterdam. However, the intervals are rather wide.

library(forecast)

par(mfrow=c(3,4),mar=c(1, 2.5, 2.5, .1) + 0.1)

sapply(unique(all$regionnum),function(x) {

      regio <- all[all$regionnum==x,]

      regio <- regio[order(regio$time),]

      top <-ts(regio$price,start=c(min(regio$year),1),

          end=c(max(regio$year),3),

          deltat=0.25)

      fit <- auto.arima(top,stepwise=FALSE,approximation=FALSE)

      LH.pred<-forecast(fit,n.ahead=4)

      plot(LH.pred,xlim=c(2000,2015),

          ylim=c(min(min(LH.pred$lower),min(regio$price[regio$year>2000])),

              max(max(LH.pred$upper),max(regio$price[regio$year>2000]))),

          main=regio$region[1])

    })

Conclusion

The prices indicate that we may indeed be close to the bottom. We are getting close to the longer term trend. However, just as they might swing over the long term trend, so they might swing under. We do see some prices bottoming out. Unfortunately, by the time we know for sure, the question if Q3 2013 is the point where it stopped getting worse may be a year or more away, at which point it is not so much of practical importance.

Appendix: reading data

vect <- function(x)   unlist(strsplit(x,' '))
dd <- dir(pattern='nvm.*pdf')
la <- lapply(dd[1:9], function(x) {
      #http://www.foolabs.com/xpdf/
      system(paste('pdftotext',x))
      lin1 <- readLines(  gsub('pdf','txt',x))
      region <- lin1[3]
      blok <- grep('^jaar',lin1,value=TRUE)
      jaar <- substr(blok,regexpr('jaar',blok)+5,regexpr('zuiver',blok)-2)
      blok <- substring(blok,regexpr('zuiver kwartaal',blok)+16)
      zk1 <- substr(blok,1,regexpr('voortschrijdend',blok)-2)
      blok <- substring(blok,regexpr('kwartaal',blok)+9)
      vk1 <- substr(blok,1,regexpr('jaar',blok)-2)
      zk2 <- substring(blok,regexpr('kwartaal',blok)+9)
      vk2 <- grep('^voortschrijdend kwartaal',lin1,value=TRUE)[1]
      vk2 <- substring(vk2,26,regexpr('jaar',vk2)-2)
      
      part1 <- data.frame(jaar=vect(jaar),
          vk1=vect(vk1),
          zk1=vect(zk1),
          vk2=vect(vk2),
          zk2=vect(zk2))
      
      blok <- grep('^periode',lin1,value=TRUE)
      jaar <- substr(blok,9,regexpr('zuiver',blok)-2)
      blok <- substring(blok,regexpr('zuiver kwartaal',blok)+16)
      zk1 <- substr(blok,1,regexpr('voortschrijdend',blok)-2)
      blok <- substring(blok,regexpr('kwartaal',blok)+9)
      vk1 <- substr(blok,1,regexpr('jaar',blok)-2)
      zk2 <- substring(blok,regexpr('kwartaal',blok)+9)
      vk2 <- grep('^voortschrijdend kwartaal',lin1,value=TRUE)[2]
      vk2 <- substring(vk2,26,regexpr('jaar',vk2)-2)
      part2 <- data.frame(jaar=vect(jaar),
          vk1=vect(vk1),
          zk1=vect(zk1),
          vk2=vect(vk2),
          zk2=vect(zk2))
      parts <- rbind(part1,part2)
      parts$region <- region
      parts
    })

# extracting data for Netherlands from region data
nl <- data.frame(
    jaar=la[[1]]$jaar,
    zk1=la[[1]]$zk2,
    region='Regio 77 Netherlands'
)
# putting together phase 1 

la2 <- lapply(la,function(x) x[,c(1,3,6)]) 
la2[[10]] <- nl
all <- do.call(rbind,la2)

# current data
new <- read.table(textConnection('
"Regio 12 Zuidwest-Drenthe" 184 189 194
"Regio 31 Waterland" 190 209 201
"Regio 34 Amsterdam" 219 224 228
"Regio 37 Almere" 168 171 161
"Regio 46 Den Haag" 213 215 202
"Regio 49 Rotterdam" 172 181 169
"Regio 55 Ede eo" 214 206 224
"Regio 65 Zeeuwse Eilanden" 184 182 184
"Regio 71 Eindhoven eo" 204 211 210
"Regio 77 Netherlands" 205 207 205
'), col.names=c('region','13-01','13-02','13-03')
)

new <- reshape(new,
    varying=list(names(new[-1])),
    v.names='price',
    timevar='jaar', 
    idvar='region',
    times=c('13-1','13-2','13-3'),
    direction='long',
)
#combining 2
all$price <- as.numeric(sub(',','.',as.character(all$zk1)))
all <- subset(all,select=-zk1)
all <- rbind(all,new)
# variable coding & conversion
all$year <- as.numeric(substr(all$jaar,1,2))
all$year <- ifelse(all$year<50,all$year+2000,all$year+1900)
all$Quarter <- substr(all$jaar,4,4)
all$time <- all$year+as.numeric(all$Quarter)/4-1/8
all$region <- factor(all$region)
all$regionnum <- as.numeric(substr(all$region,7,9))
levels(all$region) <- substr(levels(all$region),10,100)

Influence Analysis for Repeated Measures Data

I am trying exercise 59.8 (page 5057) of the SAS/STAT Users Guide 12.3 in R. The interesting thing is that influence is investigated on subject level rather than individual level. The diagnostics in nlme does not do leave-subject-out, at least, not that I know of. MCMCglm hardly has any diagnostics. This does not mean no validation is possible, this is R, programming is not optional, but rather expected. Hence with a little bit of work it is possible to estimate PRESS, Cook's D and effects on fixed effects. From this it follows extensive validation is possible, provided we can extract the underlying variables from the model fit object.

Data

Data is same as exercise 59.2 (exercise in R).
r1 <- read.table(textConnection('
1 F 21.0 20.0 21.5 23.0
2 F 21.0 21.5 24.0 25.5
3 F 20.5 24.0 24.5 26.0
4 F 23.5 24.5 25.0 26.5
5 F 21.5 23.0 22.5 23.5
6 F 20.0 21.0 21.0 22.5
7 F 21.5 22.5 23.0 25.0
8 F 23.0 23.0 23.5 24.0
9 F 20.0 21.0 22.0 21.5
10 F 16.5 19.0 19.0 19.5
11 F 24.5 25.0 28.0 28.0
12 M 26.0 25.0 29.0 31.0
13 M 21.5 22.5 23.0 26.5
14 M 23.0 22.5 24.0 27.5
15 M 25.5 27.5 26.5 27.0
16 M 20.0 23.5 22.5 26.0
17 M 24.5 25.5 27.0 28.5
18 M 22.0 22.0 24.5 26.5
19 M 24.0 21.5 24.5 25.5
20 M 23.0 20.5 31.0 26.0
21 M 27.5 28.0 31.0 31.5
22 M 23.0 23.0 23.5 25.0
23 M 21.5 23.5 24.0 28.0
24 M 17.0 24.5 26.0 29.5
25 M 22.5 25.5 25.5 26.0
26 M 23.0 24.5 26.0 30.0
27 M 22.0 21.5 23.5 25.0
'),col.names=c('Person','Gender','Age8','Age10','Age12','Age14'),
colClasses=c('factor','factor',rep('numeric',4)))
rm <- reshape(r1,direction='long',
    varying=list(c('Age8','Age10','Age12','Age14')),
    timevar='Age',idvar=c('Person','Gender'),
    v.names='y',
    times=c(8,10,12,14))
rm$Gender <- relevel(rm$Gender,ref='M')
rm$fage=factor(rm$Age)
rm$Person <- factor(rm$Person,levels=format(1:27,trim=TRUE))
rm <- rm[order(rm$Person,rm$Age),]

nlme

Analysis, standard plot are not too difficult.

lSymm <- lme(y ~  Age * Gender,

    data=rm, random= list(Person =pdSymm(~ fage-1)),method='ML')

plot(lSymm, resid(., type = "p") ~ Age | Person)

Subject level influence

I have chosen to display three items, PRESS, Cook's D and effect on fixed parameters. PRESS is reasonable straightforward, the numbers more or less match. Cook's D on the other hand, is not. I took the formula of the SAS/STAT guide, which calculated it, my wording, as Mahalanobis distance of leave-subject-out fixed parameters. However, the numbers don't match. Part of that may be that I do recalculate random parameters too. Compare my figure with Output 59.8.3 top left, this seems quite similar. 

coefFulllme <- as.numeric(coef(lSymm)[1,1:4])

VCMlme <- vcov(lSymm)

lSymmLSO <- sapply(levels(rm$Person), function(x) {

      rloo <- rm[rm$Person !=x,]

      lSymm <- lme(y ~  Age * Gender,

          data=rloo, random= list(Person =pdSymm(~ fage-1)),

          method='ML')

      coef <- as.numeric(coef(lSymm)[1,1:4])

      genderF <- rm[rm$Person==x,'Gender'][1]=='F'

      pred <- coef[1]+

          c(8,10,12,14)*coef[2]+

          genderF*coef[3]+

          c(8,10,12,14)*coef[4]*genderF

      obs <- rm[rm$Person ==x,'y']

      CD=mahalanobis(coef,coefFulllme,VCMlme)/4

      c(press=sum((obs-pred)^2),cd=CD,coef)

    })

lSymmLSO <- t(lSymmLSO)

lSymmLSO[,c(1,2)]

       press          cd

1  10.323531 0.020169407

2   3.839932 0.043571109

3  10.899537 0.032741341

4  24.050823 0.045570019

5   1.690000 0.016239351

6  11.866688 0.016490548

7   1.191555 0.005405341

8   4.678438 0.027992197

9  13.488546 0.042166906

10 85.581754 0.155169369

11 68.465309 0.106829228

12 39.449016 0.022192784

13 12.920225 0.007427751

14  6.119853 0.001987303

15 26.126968 0.277545188

16 21.065821 0.018968667

17 10.050641 0.016763371

18  7.800734 0.007611981

19 15.196358 0.008490818

20 43.256052 0.395086129

21 96.107295 0.115497649

22 13.879678 0.033635096

23  4.961204 0.011979580

24 41.758892 0.905440147

25  4.643721 0.095851759

26  8.019698 0.026621476

27 19.920691 0.017942506

plot(y=lSymmLSO[,2],x=1:27,main="Cook's D",xlab='Subject',type='h')

Fixed-Effects Deletion Estimates

Having done all the pre-work, the fixed effects deletion statistics are just a plot away.They do look slightly different from PROC MIXED as the model is a bit different in the SAS/STAT Guide.

par(mfrow=c(2,2))

dummy <- sapply(1:4,function(x) {

  plot(y=lSymmLSO[,x+2],x=1:27,main=names(coef(lSymm))[x],ylab='',xlab='Subject')

  abline(h=coefFulllme[x])

})

MCMCglm

The model is easy enough to fit. The approach used in nlme is easy enough to convert. Only first 5 subject's data shown for brevity.
library(MCMCglmm)
prior1 <- list(R=list(V=diag(4),nu=.01),
    G=list(G1=list(V=diag(1),nu=.01) ))

m1 <- MCMCglmm(y ~ Age* Gender , 
    random= ~ Person ,
    rcov=~ us(fage) :Person,
    data=rm,family='gaussian',
    nitt=500000,thin=20,burnin=50000,
    verbose=FALSE
 ,   prior=prior1
)

VCMM1 <- cov(m1$Sol)
coefFullM1 <- colMeans(m1$Sol)

m1LSO <- sapply(levels(rm$Person), function(x) {

      rloo <- rm[rm$Person !=x,]

      m1 <- MCMCglmm(y ~ Age* Gender , 

          random= ~ Person ,

          rcov=~ us(fage) :Person,

          data=rloo,family='gaussian',

     nitt=500000,thin=20,burnin=50000,

          verbose=FALSE

          ,   prior=prior1

      )

      coef <- colMeans(m1$Sol)

      genderF <- rm[rm$Person==x,'Gender'][1]=='F'

      pred <- coef[1]+

          c(8,10,12,14)*coef[2]+

          genderF*coef[3]+

          c(8,10,12,14)*coef[4]*genderF

      obs <- rm[rm$Person ==x,'y']

      CD=mahalanobis(coefFullM1,coef,VCMM1)/4

      c(press=sum((obs-pred)^2),cd=CD,coef)

    })

m1LSO <- t(m1LSO)

m1LSO[1:5,1:2]

      press         cd

1 10.133019 0.01117458

2  3.835048 0.03347745

3 10.891120 0.02740768

4 24.167881 0.03618814

5  1.694508 0.01317537

For both PRESS and Cook's D it is close. For completeness the plot.Not exactly the same, but the trend and interpretations are clear enough.

plot(y=m1LSO[,2]*4,x=1:27,main="Cook's D",xlab='Subject',type='h')

Fixed-Effects Deletion Estimates

Following the template from nlme, the result is similar. 

par(mfrow=c(2,2))

dummy <- sapply(1:4,function(x) {

      plot(y=m1LSO[,x+2],x=1:27,main=names(coefFullM1)[x],ylab='',xlab='Subject')

      abline(h=coefFullM1[x])

    })

Mixed Models: Influence in Heterogeneous Variance Model

In this post I extend my knowledge of mixed models by redoing section 59.7 (page 5048) of the SAS/STAT user guide. I don't think this particular example can be run in lme4, so that leaves nlme and MCMCglmm. MCMCglmm has less ability for influence measures, so this is only touched upon.

data

To quote the SAS/STAT User guide: 'a one-way classification model with heterogeneous variances is fit. The data, (...), represent amounts of different types of fat absorbed by batches of doughnuts during cooking, measured in grams'. Data are read below. The final statement is just ordering the data similar to SAS, which is convenient for comparison on my side. 

r1 <- read.table(textConnection('

1 164 1 172 1 168 1 177 1 156 1 195

2 178 2 191 2 197 2 182 2 185 2 177

3 175 3 193 3 178 3 171 3 163 3 176

4 155 4 166 4 149 4 164 4 170 4 168

') ,col.names=paste('a',1:12,sep=''))

r2 <- data.frame(FatType=factor(as.matrix(r1[,seq(1,11,by=2)])),

    Absorbed=as.numeric(as.matrix(r1[,seq(2,12,by=2)])))

r2 <- r2[c(seq(1,24,by=4),seq(2,24,by=4),seq(3,24,by=4),seq(4,24,by=4)),]

nlme

For a change, the gls() function is used from the nlme package. A few preparations are needed. A variable x is used as dummy variable in the varIdent sub-statement. Variable i is used in influence diagnostic plotting. SAS contrasts speak for themselves. 

r2$x <- 1

r2$i <-1:24

options(contrasts=c( 'contr.SAS','contr.SAS'))

library(nlme)

lme1 <- gls(Absorbed ~ FatType,

    weights=varIdent(form=~x |  FatType),data=r2 )    

summary(lme1)

Generalized least squares fit by REML

  Model: Absorbed ~ FatType 

  Data: r2 

       AIC      BIC    logLik

  170.3109 178.2767 -77.15543

Variance function:

 Structure: Different standard deviations per stratum

 Formula: ~x | FatType 

 Parameter estimates:

        1         2         3         4 

1.0000000 0.5825169 0.7404827 0.6162593 

Coefficients:

            Value Std.Error  t-value p-value

(Intercept)   162  3.356586 48.26332  0.0000

FatType1       10  6.397916  1.56301  0.1337

FatType2       23  4.618803  4.97965  0.0001

FatType3       14  5.247222  2.66808  0.0148

 Correlation: 

         (Intr) FtTyp1 FtTyp2

FatType1 -0.525              

FatType2 -0.727  0.381       

FatType3 -0.640  0.336  0.465

Standardized residuals:

          Min            Q1           Med            Q3           Max 

-1.581138e+00 -6.625646e-01 -6.533961e-15  5.473172e-01  1.723923e+00 

Residual standard error: 13.34166 

Degrees of freedom: 24 total; 20 residual

As ever, nlme shows results differently from SAS. The missing bit is the variances. These can be displayed by this statement.

(lme1$sigma*

      coef(lme1$modelStruct$varStruct,

          uncons = FALSE, allCoef = TRUE))^2

        1         2         3         4 

177.99996  60.39999  97.59999  67.60003 

diagnostics

Having set up the model, leave-one-out statistics is just a sapply() away. For brevity only a few are listed. 

sa <- sapply(1:nrow(r2) ,function(x) {

      rx <- r2

      rx$Absorbed[x] <- NA

      lmex <- gls(Absorbed ~ FatType,

          weights=varIdent(form=~x |  FatType),data=rx,na.action=na.omit )    

      cc <- coef(lmex$modelStruct$varStruct,uncons = FALSE, allCoef = TRUE)

      cc <- cc[order(names(cc))]

      c(coef(lmex),

         (lmex$sigma*cc)^2)

   })

tsa <- as.data.frame( t(sa))

head(tsa)

  (Intercept) FatType1 FatType2 FatType3         1        2        3        4

1         162     11.6       23       14 203.29918 60.40017 97.59971 67.60023

2         162     10.0       23       14 222.50049 60.39999 97.59995 67.59992

3         162     10.8       23       14 217.70003 60.40000 97.59999 67.59999

4         162      9.0       23       14 214.99969 60.40013 97.59964 67.60018

5         162     13.2       23       14 145.70021 60.39986 97.60000 67.60007

6         162      5.4       23       14  63.79997 60.40006 97.59999 67.59996

Plots are more interesting. 

tsa$i <- 1:nrow(tsa)

library(lattice)

par(mfrow=c(2,2))

plot(`(Intercept)`~i,data=tsa,xlab='')

plot(FatType1~i,data=tsa,xlab='')

plot(FatType2~i,data=tsa,xlab='')

plot(FatType3~i,data=tsa,xlab='')

It is a bit funny in influences, since fat-types are displayed equivalent to intercept. In terms of influence on the FatTypes means a plot not shown in SAS/STAT:

names(tsa)[1] <- 'Intercept'

xyplot(Intercept + I(Intercept+FatType1) 

        + I(Intercept+FatType2)+ I(Intercept+FatType3) ~ i   ,

    data=tsa,

    outer=TRUE,

    strip=function(..., factor.levels)

      strip.default(  ...,factor.levels=paste('FatType',c(4,1:3))))

The next plot shows influence on the estimated variances

names(tsa)[5:8] <- paste('FT',1:4,sep='')

xyplot(FT1 + FT2 + FT3 + FT4 ~  i   ,data=tsa,outer=TRUE)

residuals

The ground work for residuals is also present. Note that the studentized residuals come out slightly different; there seems to be a factor (5/6)^2 difference. 

resdf <- data.frame(Absorbed=r2$Absorbed,

    predicted=predict(lme1),

    residual=residuals(lme1),

    std=attr(residuals(lme1),'std')

    )

resdf$lofit <- 0

resdf$lofit[r2$FatType=='1'] <- 

   rowSums(tsa[,c('Intercept','FatType1')])[r2$FatType=='1']

resdf$lofit[r2$FatType=='2'] <- 

   rowSums(tsa[,c('Intercept','FatType2')])[r2$FatType=='2']

resdf$lofit[r2$FatType=='3'] <- 

   rowSums(tsa[,c('Intercept','FatType3')])[r2$FatType=='3']

resdf$lofit[r2$FatType=='4'] <- tsa[,c('Intercept')][r2$FatType=='4']

resdf$PRESSres <- resdf$Absorbed-resdf$lofit

resdf$StudRes <- resdf$residual/resdf$std

head(resdf)

   Absorbed predicted      residual      std lofit      PRESSres       StudRes

1       164       172 -8.000000e+00 13.34166 173.6 -9.600000e+00 -5.996254e-01

5       172       172 -2.842171e-14 13.34166 172.0  5.684342e-14 -2.130297e-15

9       168       172 -4.000000e+00 13.34166 172.8 -4.800000e+00 -2.998127e-01

13      177       172  5.000000e+00 13.34166 171.0  6.000000e+00  3.747659e-01

17      156       172 -1.600000e+01 13.34166 175.2 -1.920000e+01 -1.199251e+00

21      195       172  2.300000e+01 13.34166 167.4  2.760000e+01  1.723923e+00

nlme's own tools

nlme is not without tools for diagnostics as shown below. I labelled a number of points mentions in SAS/STAT Guide. The cut-off for labeling is determined by id. The value of 0.15 coincides with a two sided 85% interval, is chosen so as to label the four points mentioned in SAS/STAT Guide. It is certainly not a limit I would choose and I have some doubt about four influential points on a total of 24 points.

plot(lme1,form= resid(., type = "p") ~ i | FatType,id=.15)

Second diagnostic plot contains normal plots. The marked points seem reasonable in line with the others.

qqnorm(lme1, ~ resid(., type = "p") | FatType,id=.15)

MCMCglmm

While the nlme sections is really long, MCMCglmm is really short. To run the model there were some tricks. An additional variable avar is nested under FatTypes. MCMCglmm gives a units variance, which does not make sense in the context of this model. I did not know how to remove it, but a suitable prior was used to suppress it. 

library(MCMCglmm)

r2$avar <- factor(rep(1:6,4))

prior1 <- list(R=list(V=diag(1)*.1,nu=10),

    G=list(G1=list(V=diag(4),nu=.01) )

)

mc1 <- MCMCglmm(Absorbed ~ FatType,

    random=~ idh(FatType):avar,

    data=r2 ,pr=TRUE,verbose=FALSE,

    prior=prior1,    nitt=50000)

summary(mc1)

 Iterations = 3001:49991

 Thinning interval  = 10

 Sample size  = 4700 

 DIC: 36.68339 

 G-structure:  ~idh(FatType):avar

       post.mean l-95% CI u-95% CI eff.samp

1.avar     295.5    42.52    773.6     4700

2.avar     100.1    14.96    259.5     4712

3.avar     164.9    21.59    443.2     4700

4.avar     110.4    18.25    291.0     4700

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp

units    0.1256  0.03751   0.2588     4359

 Location effects: Absorbed ~ FatType 

            post.mean l-95% CI u-95% CI eff.samp   pMCMC    

(Intercept)  161.9363 153.0445 169.8837     4700 < 2e-04 ***

FatType1      10.0627  -6.9373  25.8132     4969 0.18681    

FatType2      22.9649  12.0911  34.9963     4700 0.00128 ** 

FatType3      13.9831   0.8961  27.4354     4700 0.04043 *  

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As before, the variances are larger than expected. There is also a function which looks dfferently at this. These values are lower than nlme & PROC MIXED.

posterior.mode(mc1$VCV)

      1.avar       2.avar       3.avar       4.avar        units 

133.71520455  42.92154371  73.72997958  51.18475597   0.09115953 

Influence diagnostics

It is easy enough to extract parameter estimates from the simulations, so a sapply() to get leave-one-out estimates similar to nlme is obvious. However,  residuals() is not yet implemented, so that is far away. It was my choice not to do this.