[R] Lmer, heteroscedasticity and permutation, need help please

[Alan Juilland]

Hi everybody,

I'm trying to analyse a set of data with a non-normal response, 2 fixed 
effects and 1 nested random effect with strong heteroscedasticity in the 
model.

I planned to use the function lmer : lmer(resp~var1*var2 + (1|rand)) and 
then use permutations based on the t-statistic given by lmer to get 
p-values.

1/ Is it a correct way to obtain p-values for my variables ? (see below)

2/ I read somewhere that lme is more adequate when heteroscedasticity is 
strong. Do I have to use lme instead of lmer ?

3/ It is possible to fit a glm in lmer using family="...". Is it 
possible to use it in spite of hard heteroscedasticity ?

4/ A last question concerning SAS. My model appears to not converge in 
SAS, indicating a "structure" in the variance. Is it implying something 
in lmer or lme ?

Many Thanks

Here is the function for the permutation :

permut<-function(data,vardep,var1,var2,var.random,nboot=10){

# initialisation
ms.var1<-numeric(length=nboot)
ms.var2<-numeric(length=nboot)
ms.var3<-numeric(length=nboot)
var1.permut<-numeric(length=length(var1))
var2.permut<-numeric(length=length(var2))
var3.permut<-numeric(length=length(var2))
var1<-factor(var1)
var2<-factor(var2)

# model initial
model.ini<-lmer(vardep~factor(var1)*factor(var2)+(1|var.random),data=data)
vc<-vcov(model.ini)
b <- fixef(model.ini)
se<- sqrt(diag(vc))
z<- b/se
p<- 2*(1-pnorm(abs(z)))

# boucle de permutation
for(i in 1:nboot){

#boucles de randomisation des traitements
for (j in 1:nlevels(var1)){
var2.permut[var1==j]<-sample(var2[var1==j])}
for (j in 1:nlevels(var2)){
var1.permut[var2==j]<-sample(var1[var2==j])}
data2<-data.frame(vardep,var1,var2,var.random,var1.permut,var2.permut)

#test var1:
m1<-lmer(vardep~factor(var1.permut)*factor(var2)+(1|var.random),data=data2)
ms.var1[i]<-as.numeric(fixef(m1)/sqrt(diag(vcov(m1))))[2]
#test var2:
m2<-lmer(vardep~factor(var1)*factor(var2.permut)+(1|var.random),data=data2)
ms.var2[i]<-as.numeric(fixef(m2)/sqrt(diag(vcov(m2))))[3]
#test interaction:
m3<-lmer(vardep~factor(var1.permut)*factor(var2.permut)+(1|var.random),data=data2)
ms.var3[i]<-as.numeric(fixef(m3)/sqrt(diag(vcov(m3))))[4]
}

max.boot<-c("NA",max(ms.var1),max(ms.var2),max(ms.var3))
min.boot<-c("NA",min(ms.var1),min(ms.var2),min(ms.var3))
pval<-c("NA",length(ms.var1[ms.var1>=z[2]])/nboot,length(ms.var2[ms.var2>=z[3]])/nboot,length(ms.var3[ms.var3>=z[4]])/nboot)

res<-cbind(b,se,z,p,permut.min=min.boot,permut.max=max.boot,permut.pval=pval)
par(mfrow=c(2,2))
hist(ms.var1);abline(v=z[2],col="RED");hist(ms.var2);abline(v=z[3],col="RED");hist(ms.var3);abline(v=z[4],col="RED")

return(res)
}



-- 
Alan Juilland – PhD Student
Department of Ecology and Evolution
Biophore, University of Lausanne
1015 Dorigny
Switzerland
Tel : ++41 21 692 41 74
Fax :  +41 21 692 41 65