[R-sig-ME] VPC calculation

[Dr C B Stride]

Hi
Can anyone point me towards some R code, and pref an example or two of 
using it, for calculating the Variance Partition Coefficient (VPC) for a 
multilevel logistic regression and also a multilevel poisson regression 
using the simulation method outlined in Goldstein et al (2002)? cheers
Chris


John.Morrongiello at csiro.au said the following on 01/06/2011 08:17:
> Hi
> Is there a way to investigate auto correlation in a model fit with lmer, similar to the ACF.lme function for lme models? I haven't been able to find any thread discussion on the topic and some advice would be much appreciated. I thought that maybe I could use the by-subject outputs from acf.fnc  (languageR package), but this doesn't work (see below).
> 
> library(nlme)
> library(lme4)
> library(languageR)
> 
> ##consider this data set
> nID<-20
> nwithin<-10
> prosw<-nID*nwithin
> 
> growth<-vector(length=prosw)
> age<-seq(2,11, length=nwithin)
> a<-rnorm(prosw,mean=2)
> b<-120;c<-10;d<-5;e<- -2;rho<--0.9
> for (i in 1:nID)
> {
> sf<-0.1*rnorm(1);sc<-2*rnorm(1);ss<-rnorm(1,mean=2);err<-0.1*rnorm(1)
> for (j in 1:nwithin)
> {
> growth[nwithin*(i-1)+j]<-sf +  b + (c+sc)*a[nwithin*(i-1)+j] - (d+ss)*age[j] +sf*a[nwithin*(i-1)+j]+
>                    err
> err<-err*rho + 0.1*rnorm(1)
> }
> }
> a<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> b<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> c<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> d<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> e<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> f<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> g<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> h<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> i<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> j<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> k<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> l<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> m<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> n<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> o<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> p<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> q<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> r<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> s<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> t<-seq((1980+sample(0:10,1)),length.out=10,by=1)
> year<-c(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t)##sorry couldn't get for loop to work for this!
> age<-rep(age, nID)
> ID<-rep(seq(1,nID), each=nwithin)
> growthdata<-data.frame(ID,age,year,growth)
> growthdata$ID<-as.factor(growthdata$ID)
> 
> ##I would like to fit a mixed model with crossed random effects of ID and year (not every level of ID is present in each year). It is highly likely that there is correlation between sequential measurements within ID. I can account for this using an auto-regressive correlation structure of the form correlation=corAR1(form=~age|ID) in lme, but due to the high number of ID and year levels in my full dataset, crossing ID with year is not possible. Therefore, I thought I'd try lmer and fit the model:
> 
> model.lmer<-lmer(growth~age+(age|ID)+(1|as.factor(year)),growthdata)
> 
> ###  and see how if performed. I fit a similar model in lme to compare the outputs of ACF and acf.fnc using model residuals. I naively thought that averaging across levels of 'lag' in the acf.fnc output may work, but it obviously doesn't! My real data has a much higher level of autocorrelation.
> 
> model.lme<-lme(growth~age,random=~age|ID,growthdata)
> ACF(model.lme)
> lag         ACF
>    0  1.00000000
>    1 -0.16324159
>    2 -0.05552011
>    3 -0.16254222
>    4 -0.10331816
>    5 -0.06783792
>    6 -0.03268318
>    7  0.08196762
>    8 -0.34495610
>    9  0.17929713
> growthdata$resid<-resid(model.lme)
> y<-acf.fnc(growthdata,group='ID',time='age',x='resid',plot=F)
> yy<-as.matrix(tapply(y$Acf,y$Lag,mean))
>          [,1]
> 0  1.00000000
> 1 -0.12350304
> 2 -0.10730408
> 3 -0.04293020
> 4 -0.09515590
> 5 -0.07507144
> 6  0.01232782
> 7  0.01487862
> 8 -0.07217361
> 9 -0.01106817
> 
> Cheers
> 
> John
> 
> 
> 
> 
> 
> 
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


-- 
Dr Chris Stride, C. Stat, Statistician, Institute of Work Psychology, 
University of Sheffield
Telephone: 0114 2223262
Fax: 0114 2727206

“Figure It Out”
Statistical Consultancy and Training Service for Social Scientists

Visit www.figureitout.org.uk for details of my consultancy services, and 
forthcoming training courses, which are also available on an in-house basis:
- Data Management using SPSS syntax
- Multiple Regression using SPSS
- Multilevel Modelling using SPSS
- Structural Equation Modelling using MPlus