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
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
