Lorsque vous faites des inférences sur les moyennes de groupe, les intervalles crédibles sont-ils sensibles à la variance intra-sujet alors que les intervalles de confiance ne le sont pas?

8

Ceci est un spin-off de cette question: comment comparer deux groupes avec plusieurs mesures pour chaque individu avec R?

Dans les réponses (si j'ai bien compris), j'ai appris que la variance intra-sujet n'affecte pas les inférences faites sur les moyennes de groupe et il est correct de simplement prendre les moyennes des moyennes pour calculer la moyenne de groupe, puis de calculer la variance intra-groupe et de l'utiliser pour effectuer des tests de signification. Je voudrais utiliser une méthode où plus la variance intra-sujet est grande, moins je suis sûr du groupe signifie ou comprendre pourquoi il n'est pas logique de désirer cela.

Voici un tracé des données originales avec quelques données simulées qui ont utilisé les mêmes moyennes de sujet, mais ont échantillonné les mesures individuelles pour chaque sujet à partir d'une distribution normale en utilisant ces moyennes et une petite variance intra-sujet (sd = 0,1). Comme on peut le voir, les intervalles de confiance au niveau du groupe (ligne du bas) ne sont pas affectés par cela (du moins la façon dont je les ai calculés).

entrez la description de l'image ici

J'ai également utilisé rjags pour estimer les moyennes du groupe de trois manières. 1) Utiliser les données originales brutes 2) Utiliser uniquement les moyens Sujet 3) Utiliser les données simulées avec un petit sd intra-sujet

Les résultats sont ci-dessous. En utilisant cette méthode, nous voyons que les intervalles crédibles à 95% sont plus étroits dans les cas n ° 2 et n ° 3. Cela correspond à mon intuition de ce que j'aimerais me produire lorsque je fais des inférences sur les moyens de groupe, mais je ne suis pas sûr que ce soit juste un artefact de mon modèle ou une propriété d'intervalles crédibles.

Remarque. Pour utiliser rjags, vous devez d'abord installer JAGS à partir d'ici: http://sourceforge.net/projects/mcmc-jags/files/

entrez la description de l'image ici

Les différents codes sont ci-dessous.

Les données d'origine:

structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 
3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 
6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 
10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 
12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15, 
15, 15, 15, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 18, 
18, 18, 18, 18, 18, 2, 0, 16, 2, 16, 2, 8, 10, 8, 6, 4, 4, 8, 
22, 12, 24, 16, 8, 24, 22, 6, 10, 10, 14, 8, 18, 8, 14, 8, 20, 
6, 16, 6, 6, 16, 4, 2, 14, 12, 10, 4, 10, 10, 8, 4, 10, 16, 16, 
2, 8, 4, 0, 0, 2, 16, 10, 16, 12, 14, 12, 8, 10, 12, 8, 14, 8, 
12, 20, 8, 14, 2, 4, 8, 16, 10, 14, 8, 14, 12, 8, 14, 4, 8, 8, 
10, 4, 8, 20, 8, 12, 12, 22, 14, 12, 26, 32, 22, 10, 16, 26, 
20, 12, 16, 20, 18, 8, 10, 26), .Dim = c(108L, 3L), .Dimnames = list(
    NULL, c("Group", "Subject", "Value")))

Obtenir des moyennes de sujet et simuler les données avec une petite variance intra-sujet:

#Get Subject Means
means<-aggregate(Value~Group+Subject, data=dat, FUN=mean)

#Initialize "dat2" dataframe
dat2<-dat

#Sample individual measurements for each subject
temp=NULL
for(i in 1:nrow(means)){
  temp<-c(temp,rnorm(6,means[i,3], .1))
}

#Set Simulated values
dat2[,3]<-temp

La fonction pour s'adapter au modèle JAGS:

 require(rjags) 

#Jags fit function
jags.fit<-function(dat2){

  #Create JAGS model
  modelstring = "

  model{
  for(n in 1:Ndata){
  y[n]~dnorm(mu[subj[n]],tau[subj[n]]) T(0, )
  }

  for(s in 1:Nsubj){
  mu[s]~dnorm(muG,tauG) T(0, )
  tau[s] ~ dgamma(5,5)
  }


  muG~dnorm(10,.01) T(0, )
  tauG~dgamma(1,1)

  }
  "
  writeLines(modelstring,con="model.txt")

#############  

  #Format Data
  Ndata = nrow(dat2)
  subj = as.integer( factor( dat2$Subject ,
                             levels=unique(dat2$Subject ) ) )
  Nsubj = length(unique(subj))
  y = as.numeric(dat2$Value)

  dataList = list(
    Ndata = Ndata ,
    Nsubj = Nsubj ,
    subj = subj ,
    y = y
  )

  #Nodes to monitor
  parameters=c("muG","tauG","mu","tau")


  #MCMC Settings
  adaptSteps = 1000             
  burnInSteps = 1000            
  nChains = 1                   
  numSavedSteps= nChains*10000          
  thinSteps=20                      
  nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains )            


  #Create Model
  jagsModel = jags.model( "model.txt" , data=dataList, 
                          n.chains=nChains , n.adapt=adaptSteps , quiet=FALSE )
  # Burn-in:
  cat( "Burning in the MCMC chain...\n" )
  update( jagsModel , n.iter=burnInSteps )

  # Getting DIC data:
  load.module("dic")


  # The saved MCMC chain:
  cat( "Sampling final MCMC chain...\n" )
  codaSamples = coda.samples( jagsModel , variable.names=parameters , 
                              n.iter=nPerChain , thin=thinSteps )  

  mcmcChain = as.matrix( codaSamples )

  result = list(codaSamples=codaSamples, mcmcChain=mcmcChain)

}

Ajustez le modèle à chaque groupe de chaque jeu de données:

#Fit to raw data
groupA<-jags.fit(dat[which(dat[,1]==1),])
groupB<-jags.fit(dat[which(dat[,1]==2),])
groupC<-jags.fit(dat[which(dat[,1]==3),])

#Fit to subject mean data
groupA2<-jags.fit(means[which(means[,1]==1),])
groupB2<-jags.fit(means[which(means[,1]==2),])
groupC2<-jags.fit(means[which(means[,1]==3),])

#Fit to simulated raw data (within-subject sd=.1)
groupA3<-jags.fit(dat2[which(dat2[,1]==1),])
groupB3<-jags.fit(dat2[which(dat2[,1]==2),])
groupC3<-jags.fit(dat2[which(dat2[,1]==3),])

Fonction d'intervalle crédible / intervalle de densité la plus élevée:

#HDI Function
get.HDI<-function(sampleVec,credMass){ 
  sortedPts = sort( sampleVec )
  ciIdxInc = floor( credMass * length( sortedPts ) )
  nCIs = length( sortedPts ) - ciIdxInc
  ciWidth = rep( 0 , nCIs )
  for ( i in 1:nCIs ) {
    ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
  }
  HDImin = sortedPts[ which.min( ciWidth ) ]
  HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
  HDIlim = c( HDImin , HDImax, credMass )
  return( HDIlim )
}

Premier tracé:

layout(matrix(c(1,1,2,2,3,4),nrow=3,ncol=2, byrow=T))

boxplot(dat[,3]~dat[,2], 
xlab="Subject", ylab="Value", ylim=c(0, 1.2*max(dat[,3])),
col=c(rep("Red",length(which(dat[,1]==unique(dat[,1])[1]))/6),
rep("Green",length(which(dat[,1]==unique(dat[,1])[2]))/6),
rep("Blue",length(which(dat[,1]==unique(dat[,1])[3]))/6)
),
main="Original Data"
)
stripchart(dat[,3]~dat[,2], vert=T, add=T, pch=16)
legend("topleft", legend=c("Group A", "Group B", "Group C", "Individual Means +/- 95% CI"),
col=c("Red","Green","Blue", "Grey"), lwd=3, bty="n", pch=c(15),
pt.cex=c(rep(0.1,3),1),
ncol=3)

for(i in 1:length(unique(dat[,2]))){
  m<-mean(examp[which(dat[,2]==unique(dat[,2])[i]),3])
  ci<-t.test(dat[which(dat[,2]==unique(dat[,2])[i]),3])$conf.int[1:2]

  points(i-.3,m, pch=15,cex=1.5, col="Grey")
  segments(i-.3,
           ci[1],i-.3,
           ci[2], lwd=4, col="Grey"
  )
}



boxplot(dat2[,3]~dat2[,2], 
xlab="Subject", ylab="Value", ylim=c(0, 1.2*max(dat2[,3])),
col=c(rep("Red",length(which(dat2[,1]==unique(dat2[,1])[1]))/6),
rep("Green",length(which(dat2[,1]==unique(dat2[,1])[2]))/6),
rep("Blue",length(which(dat2[,1]==unique(dat2[,1])[3]))/6)
),
main=c("Simulated Data", "Same Subject Means but Within-Subject SD=.1")
)
stripchart(dat2[,3]~dat2[,2], vert=T, add=T, pch=16)
legend("topleft", legend=c("Group A", "Group B", "Group C", "Individual Means +/- 95% CI"),
col=c("Red","Green","Blue", "Grey"), lwd=3, bty="n", pch=c(15),
pt.cex=c(rep(0.1,3),1),
ncol=3)

for(i in 1:length(unique(dat2[,2]))){
  m<-mean(examp[which(dat2[,2]==unique(dat2[,2])[i]),3])
  ci<-t.test(dat2[which(dat2[,2]==unique(dat2[,2])[i]),3])$conf.int[1:2]

  points(i-.3,m, pch=15,cex=1.5, col="Grey")
  segments(i-.3,
           ci[1],i-.3,
           ci[2], lwd=4, col="Grey"
  )
}


means<-aggregate(Value~Group+Subject, data=dat, FUN=mean)

boxplot(means[,3]~means[,1], col=c("Red","Green","Blue"),
ylim=c(0,1.2*max(means[,3])), ylab="Value", xlab="Group",
main="Original Data"
)
stripchart(means[,3]~means[,1], pch=16, vert=T, add=T)

for(i in 1:length(unique(means[,1]))){
  m<-mean(means[which(means[,1]==unique(means[,1])[i]),3])
  ci<-t.test(means[which(means[,1]==unique(means[,1])[i]),3])$conf.int[1:2]

  points(i-.3,m, pch=15,cex=1.5, col="Grey")
  segments(i-.3,
           ci[1],i-.3,
           ci[2], lwd=4, col="Grey"
  )
}
legend("topleft", legend=c("Group Means +/- 95% CI"), bty="n", pch=15, lwd=3, col="Grey")


means2<-aggregate(Value~Group+Subject, data=dat2, FUN=mean)

boxplot(means2[,3]~means2[,1], col=c("Red","Green","Blue"),
ylim=c(0,1.2*max(means2[,3])), ylab="Value", xlab="Group",
main="Simulated Data Group Averages"
)
stripchart(means2[,3]~means2[,1], pch=16, vert=T, add=T)

for(i in 1:length(unique(means2[,1]))){
  m<-mean(means[which(means2[,1]==unique(means2[,1])[i]),3])
  ci<-t.test(means[which(means2[,1]==unique(means2[,1])[i]),3])$conf.int[1:2]

  points(i-.3,m, pch=15,cex=1.5, col="Grey")
  segments(i-.3,
           ci[1],i-.3,
           ci[2], lwd=4, col="Grey"
  )
}
legend("topleft", legend=c("Group Means +/- 95% CI"), bty="n", pch=15, lwd=3,   col="Grey")

Deuxième tracé:

layout(matrix(c(1,2,3,4,4,4,5,5,5,6,6,6),nrow=4,ncol=3, byrow=T))

#Plot priors
plot(seq(0,10,by=.01),dgamma(seq(0,10,by=.01),5,5), type="l", lwd=4,
     xlab="Value", ylab="Density",
     main="Prior on Within-Subject Precision"
)
plot(seq(0,10,by=.01),dgamma(seq(0,10,by=.01),1,1), type="l", lwd=4,
     xlab="Value", ylab="Density",
     main="Prior on Within-Group Precision"
)
plot(seq(0,300,by=.01),dnorm(seq(0,300,by=.01),10,100), type="l", lwd=4,
     xlab="Value", ylab="Density",
     main="Prior on Group Means"
)


#Set overall xmax value
x.max<-1.1*max(groupA$mcmcChain[,"muG"],groupB$mcmcChain[,"muG"],groupC$mcmcChain[,"muG"],
               groupA2$mcmcChain[,"muG"],groupB2$mcmcChain[,"muG"],groupC2$mcmcChain[,"muG"],
               groupA3$mcmcChain[,"muG"],groupB3$mcmcChain[,"muG"],groupC3$mcmcChain[,"muG"]
)


#Plot result for raw data
#Set ymax
y.max<-1.1*max(density(groupA$mcmcChain[,"muG"])$y,density(groupB$mcmcChain[,"muG"])$y,density(groupC$mcmcChain[,"muG"])$y)

plot(density(groupA$mcmcChain[,"muG"]),xlim=c(0,x.max), 
     ylim=c(-.1*y.max,y.max), lwd=3, col="Red",
     main="Group Mean Estimates: Fit to Raw Data", xlab="Value"
)
lines(density(groupB$mcmcChain[,"muG"]), lwd=3, col="Green")
lines(density(groupC$mcmcChain[,"muG"]), lwd=3, col="Blue")

hdi<-get.HDI(groupA$mcmcChain[,"muG"], .95)
segments(hdi[1],-.033*y.max,hdi[2],-.033*y.max, lwd=3, col="Red")

hdi<-get.HDI(groupB$mcmcChain[,"muG"], .95)
segments(hdi[1],-.066*y.max,hdi[2],-.066*y.max, lwd=3, col="Green")

hdi<-get.HDI(groupC$mcmcChain[,"muG"], .95)
segments(hdi[1],-.099*y.max,hdi[2],-.099*y.max, lwd=3, col="Blue")

####

#Plot result for mean data

#x.max<-1.1*max(groupA2$mcmcChain[,"muG"],groupB2$mcmcChain[,"muG"],groupC2$mcmcChain[,"muG"])
y.max<-1.1*max(density(groupA2$mcmcChain[,"muG"])$y,density(groupB2$mcmcChain[,"muG"])$y,density(groupC2$mcmcChain[,"muG"])$y)

plot(density(groupA2$mcmcChain[,"muG"]),xlim=c(0,x.max), 
     ylim=c(-.1*y.max,y.max), lwd=3, col="Red",
     main="Group Mean Estimates: Fit to Subject Means", xlab="Value"
)
lines(density(groupB2$mcmcChain[,"muG"]), lwd=3, col="Green")
lines(density(groupC2$mcmcChain[,"muG"]), lwd=3, col="Blue")

hdi<-get.HDI(groupA2$mcmcChain[,"muG"], .95)
segments(hdi[1],-.033*y.max,hdi[2],-.033*y.max, lwd=3, col="Red")

hdi<-get.HDI(groupB2$mcmcChain[,"muG"], .95)
segments(hdi[1],-.066*y.max,hdi[2],-.066*y.max, lwd=3, col="Green")

hdi<-get.HDI(groupC2$mcmcChain[,"muG"], .95)
segments(hdi[1],-.099*y.max,hdi[2],-.099*y.max, lwd=3, col="Blue")




####
#Plot result for simulated data
#Set ymax
#x.max<-1.1*max(groupA3$mcmcChain[,"muG"],groupB3$mcmcChain[,"muG"],groupC3$mcmcChain[,"muG"])
y.max<-1.1*max(density(groupA3$mcmcChain[,"muG"])$y,density(groupB3$mcmcChain[,"muG"])$y,density(groupC3$mcmcChain[,"muG"])$y)

plot(density(groupA3$mcmcChain[,"muG"]),xlim=c(0,x.max), 
     ylim=c(-.1*y.max,y.max), lwd=3, col="Red",
     main=c("Group Mean Estimates: Fit to Simulated data", "(Within-Subject SD=0.1)"), xlab="Value"
)
lines(density(groupB3$mcmcChain[,"muG"]), lwd=3, col="Green")
lines(density(groupC3$mcmcChain[,"muG"]), lwd=3, col="Blue")

hdi<-get.HDI(groupA3$mcmcChain[,"muG"], .95)
segments(hdi[1],-.033*y.max,hdi[2],-.033*y.max, lwd=3, col="Red")

hdi<-get.HDI(groupB3$mcmcChain[,"muG"], .95)
segments(hdi[1],-.066*y.max,hdi[2],-.066*y.max, lwd=3, col="Green")

hdi<-get.HDI(groupC3$mcmcChain[,"muG"], .95)
segments(hdi[1],-.099*y.max,hdi[2],-.099*y.max, lwd=3, col="Blue")

EDITER avec ma version personnelle de la réponse de @ StéphaneLaurent

J'ai utilisé le modèle qu'il a décrit pour échantillonner à partir d'une distribution normale avec une moyenne = 0, entre la variance du sujet = 1 et l'erreur / variance au sein du sujet = 0,1,1,10,100. Un sous-ensemble des intervalles de confiance est indiqué dans les panneaux de gauche tandis que la distribution de leurs largeurs est indiquée par les panneaux de droite correspondants. Cela m'a convaincu qu'il a 100% raison. Cependant, je suis toujours confus par mon exemple ci-dessus, mais je poursuivrai avec une nouvelle question plus ciblée.

entrez la description de l'image ici

Le code pour la simulation et les graphiques ci-dessus:

dev.new()
par(mfrow=c(4,2))


num.sims<-10000
sigmaWvals<-c(.1,1,10,100)
muG<-0  #Grand Mean
sigma.between<-1  #Between Experiment sd

for(sigma.w in sigmaWvals){

  sigma.within<-sigma.w #Within Experiment sd

  out=matrix(nrow=num.sims,ncol=2)
  for(i in 1:num.sims){

    #Sample the three experiment means (mui, i=1:3)
    mui<-rnorm(3,muG,sigma.between)

    #Sample the three obersvations for each experiment (muij, i=1:3, j=1:3)
    y1j<-rnorm(3,mui[1],sigma.within)
    y2j<-rnorm(3,mui[2],sigma.within)
    y3j<-rnorm(3,mui[3],sigma.within)


    #Put results in data frame
    d<-as.data.frame(cbind(
      c(rep(1,3),rep(2,3),rep(3,3)),
      c(y1j, y2j, y3j )
    ))
    d[,1]<-as.factor(d[,1])

    #Calculate means for each experiment
    dmean<-aggregate(d[,2]~d[,1], data=d, FUN=mean)

    #Add new confidence interval data to output
    out[i,]<-t.test(dmean[,2])$conf.int[1:2]

  }

  #Calculate % of intervals that contained muG
  cover<-matrix(nrow=nrow(out),ncol=1)
  for(i in 1:nrow(out)){
    cover[i]<-out[i,1]<muG & out[i,2]>muG
  }



  sub<-floor(seq(1,nrow(out),length=100))
  plot(out[sub,1], ylim=c(min(out[sub,1]),max(out[sub,2])),
       xlab="Simulation #", ylab="Value", xaxt="n",
       main=c(paste("# of Sims=",num.sims),
              paste("% CIs Including muG=",100*round(length(which(cover==T))/nrow(cover),3)))
  )
  axis(side=1, at=1:100, labels=sub)
  points(out[sub,2])

  cnt<-1
  for(i in sub){
    segments(cnt, out[i,1],cnt,out[i,2])
    cnt<-cnt+1
  }
  abline(h=0, col="Red", lwd=3)

  hist(out[,2]-out[,1], freq=F, xlab="Width of 95% CI",
       main=c(paste("muG=", muG), 
              paste("Sigma Between=",sigma.between), 
              paste("Sigma Within=",sigma.within))
  )

}
r confidence-interval mixed-model jags error-propagation Ballon
la source
Eh bien, je viens de trouver cette question. Aucune réponse n'a été acceptée: stats.stackexchange.com/questions/12002/…
Flask
C'est curieux que personne ici ne semble connaître mon "truc". Je viens de répondre à cette question.
Stéphane Laurent
Je viens de jeter un coup d'œil à votre modèle JAGS. Il est différent du modèle fréquentiste car vous supposez une variance différente pour chaque sujet (imbriqué dans le groupe).
Stéphane Laurent
... et votre modèle JAGS suppose également une différence entre les différents pour chaque groupe (parce que vous exécutez le modèle séparément pour chaque groupe, si je comprends bien)
Stéphane Laurent
1
L'astuce consiste à réduire le modèle mixte à un modèle simple en prenant pour observations les moyens des sujets dans votre cas et les moyens des groupes dans l'autre question. Je ne sais pas ce que vous devez faire mais j'ai affirmé que la distribution d'échantillonnage de votre modèle bayésien n'est pas la même que celle du modèle fréquentiste.
Stéphane Laurent

Réponses:

4

Dans les réponses (si j'ai bien compris), j'ai appris que la variance intra-sujet n'affecte pas les inférences faites sur les moyennes de groupe et il est correct de simplement prendre les moyennes des moyennes pour calculer la moyenne de groupe, puis de calculer la variance intra-groupe et de l'utiliser pour effectuer des tests de signification.

Permettez-moi de développer cette idée ici. Le modèle des observations individuelles est

yijk=μi+αij+ϵijk
, où :

  • yijk est le k-ème mesure de l'individu j du groupe i

  • αij∼iidN(0,σb2) est l'effet aléatoire pour l'individu j du groupe i

  • ϵijk∼iidN(0,σw2) est l'erreur intra-erreur

Dans ma réponse à votre première question , je vous propose de noter que l'on obtient un modèle linéaire gaussien classique (effets fixes) pour les sujets moyensy¯ij∙. En effet, vous pouvez facilement vérifier que

y¯ij∙=μi+δij
avec
δij=αij+1K∑kϵijk∼iidN(0,σ2)where σ2=σb2+σw2K,
en supposant Kmesures répétées pour chaque individu. Ce n'est rien d'autre que le modèle ANOVA unidirectionnel avec un facteur fixe.

Et puis j'ai affirmé que pour tirer des conclusions sur la μi vous pouvez simplement considérer le modèle linéaire classique simple dont les observations sont les moyens des sujets y¯ij∙. Mise à jour 12/04/2014 : Quelques exemples de cette idée sont maintenant écrits sur mon blog: Réduire un modèle pour obtenir des intervalles de confiance . J'ai l'impression que cela fonctionne toujours lorsque nous faisons la moyenne des données sur les niveaux d'un effet aléatoire.

Dans les réponses (si j'ai bien compris), j'ai appris que la variance intra-sujet n'affecte pas les inférences faites sur les moyennes de groupe et il est correct de simplement prendre les moyennes des moyennes pour calculer la moyenne de groupe, puis de calculer la variance intra-groupe et de l'utiliser pour effectuer des tests de signification. Je voudrais utiliser une méthode où plus la variance intra-sujet est grande, moins je suis sûr du groupe signifie ou comprendre pourquoi il n'est pas logique de désirer cela.

Comme vous le voyez dans la formule encadrée, la variance intra σw2 joue un rôle dans le modèle des moyennes de groupe observées.

Stéphane Laurent
la source
Ce modèle est logique, mais le code ne semble pas incorporer ces informations et je suppose que c'est ce que j'aimerais comprendre. Suis-je en train de faire cela correctement pour l'estimation des paramètres?: lmer(Value~Group -1 + (1|Subject), dat) lmer(Value~Group -1 + (1|Subject), dat2)Où dat est les données d'origine et dat2 est la simulation avec une petite variance intra-sujet. J'obtiens les mêmes erreurs standard.
Flask
Je n'ai pas essayé, mais cela semble étrange, vous supprimez l'interception fixe mais il y a une interception aléatoire par sujet. D'un point de vue théorique, je ne vois aucun problème, mais je ne sais pas exactement comment lmertraite les modèles sans intercept. Gardez l'interception pour être sûr.
Stéphane Laurent
J'ai suivi cette instruction car je ne pouvais pas autrement trouver comment obtenir une estimation d'intervalle. Ma compréhension de la syntaxe de la formule R est faible, alors cela n'a peut-être aucun sens.
Flask
@Flask AFAIK, il n'y a actuellement aucun package dans R fournissant un moyen d'obtenir des intervalles de confiance "corrects" pour les lmermodèles. Pour votre modèle dans le cas particulier d'une conception équilibrée, il existe des méthodes exactes des moindres carrés, mais je ne sais pas si elles sont disponibles dans certains packages.
Stéphane Laurent
1
Bien que je me demande si le lsmeanspackage avec le pbkrtestpackage pourrait fournir de bons intervalles de confiance.
Stéphane Laurent