Sind glaubwürdige Intervalle bei Rückschlüssen auf Gruppenmittelwerte empfindlich gegenüber Abweichungen innerhalb des Subjekts, während dies bei Konfidenzintervallen nicht der Fall ist?

8

Dies ist ein Nebeneffekt dieser Frage: Wie kann man zwei Gruppen mit mehreren Messungen für jedes Individuum mit R vergleichen?

In den Antworten dort (wenn ich richtig verstanden habe) habe ich gelernt, dass die Varianz innerhalb des Subjekts keine Rückschlüsse auf Gruppenmittelwerte hat, und es ist in Ordnung, einfach die Durchschnittswerte der Mittelwerte zur Berechnung des Gruppenmittelwerts zu verwenden, dann die Varianz innerhalb der Gruppe zu berechnen und diese zu verwenden Signifikanztests durchführen. Ich möchte eine Methode verwenden, bei der je größer die Varianz innerhalb des Subjekts ist, desto weniger sicher bin ich, was die Gruppenmittel angeht oder warum es keinen Sinn macht, dies zu wünschen.

Hier ist eine grafische Darstellung der Originaldaten zusammen mit einigen simulierten Daten, die das gleiche Subjektmittel verwendeten, aber die einzelnen Messungen für jedes Subjekt aus einer Normalverteilung unter Verwendung dieser Mittelwerte und einer kleinen Varianz innerhalb des Subjekts (sd = 0,1) abtasteten. Wie zu sehen ist, bleiben die Konfidenzintervalle auf Gruppenebene (untere Reihe) davon unberührt (zumindest so, wie ich sie berechnet habe).

Geben Sie hier die Bildbeschreibung ein

Ich habe auch rjags verwendet, um die Gruppenmittelwerte auf drei Arten zu schätzen. 1) Verwenden Sie die rohen Originaldaten. 2) Verwenden Sie nur die Subjektmittel. 3) Verwenden Sie die simulierten Daten mit kleinen sd innerhalb des Subjekts

Die Ergebnisse sind unten. Mit dieser Methode sehen wir, dass die zu 95% glaubwürdigen Intervalle in den Fällen 2 und 3 enger sind. Dies entspricht meiner Vorstellung davon, was ich gerne tun würde, wenn ich Rückschlüsse auf Gruppenmittel ziehen würde, aber ich bin mir nicht sicher, ob dies nur ein Artefakt meines Modells oder eine Eigenschaft glaubwürdiger Intervalle ist.

Hinweis. Um rjags verwenden zu können, müssen Sie zuerst JAGS von hier aus installieren: http://sourceforge.net/projects/mcmc-jags/files/

Geben Sie hier die Bildbeschreibung ein

Die verschiedenen Codes sind unten.

Die Originaldaten:

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")))

Erhalten Sie Subjektmittel und simulieren Sie die Daten mit geringer Varianz innerhalb des Subjekts:

#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

Die Funktion für das JAGS-Modell:

 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)

}

Passen Sie das Modell an jede Gruppe jedes Datensatzes an:

#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),])

Glaubwürdiges Intervall / Intervallfunktion mit der höchsten Dichte:

#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 )
}

Erste Handlung:

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")

Zweite Handlung:

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")

BEARBEITEN Sie mit meiner persönlichen Version der Antwort von @ StéphaneLaurent

Ich habe das von ihm beschriebene Modell verwendet, um aus einer Normalverteilung mit Mittelwert = 0, zwischen Subjektvarianz = 1 und innerhalb des Subjektfehlers / Varianz = 0,1,1,10,100 eine Stichprobe zu erstellen. Eine Teilmenge der Konfidenzintervalle wird in den linken Feldern angezeigt, während die Verteilung ihrer Breiten in den entsprechenden rechten Feldern angezeigt wird. Das hat mich überzeugt, dass er 100% richtig ist. Ich bin jedoch immer noch verwirrt von meinem obigen Beispiel, werde aber daraufhin eine neue, fokussiertere Frage stellen.

Geben Sie hier die Bildbeschreibung ein

Der Code für die obige Simulation und Diagramme:

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))
  )

}
Flasche
quelle
Nun, ich habe gerade diese Frage gefunden. Es wurde keine Antwort akzeptiert: stats.stackexchange.com/questions/12002/…
Flask
Es ist merkwürdig, dass hier niemand meinen "Trick" zu kennen scheint. Ich habe diese Frage gerade beantwortet.
Stéphane Laurent
Ich habe mir gerade Ihr JAGS-Modell angesehen. Es unterscheidet sich vom Frequentist-Modell, da Sie für jedes Thema (in Gruppe verschachtelt) eine andere Varianz annehmen.
Stéphane Laurent
... und Ihr JAGS-Modell geht auch von einer unterschiedlichen Zwischenvarianz für jede Gruppe aus (da Sie das Modell meines Wissens für jede Gruppe separat ausführen)
Stéphane Laurent
1
Der "Trick" besteht darin, das gemischte Modell auf ein einfaches Modell zu reduzieren, indem für Beobachtungen die Subjektmittel in Ihrem Fall und die Gruppenmittel in der anderen Frage herangezogen werden. Ich weiß nicht, was Sie tun sollen, aber ich habe behauptet, dass die Stichprobenverteilung Ihres Bayes'schen Modells nicht mit der des frequentistischen Modells übereinstimmt.
Stéphane Laurent

Antworten:

4

In den Antworten dort (wenn ich richtig verstanden habe) habe ich gelernt, dass die Varianz innerhalb des Subjekts keine Rückschlüsse auf Gruppenmittelwerte hat, und es ist in Ordnung, einfach die Durchschnittswerte der Mittelwerte zur Berechnung des Gruppenmittelwerts zu verwenden, dann die Varianz innerhalb der Gruppe zu berechnen und diese zu verwenden Signifikanztests durchführen.

yijk=μi+αij+ϵijk
  • yijkkji

  • αijiidN(0,σb2)ji

  • ϵijkiidN(0,σw2)

y¯ij

y¯ij=μi+δij
δij=αij+1KkϵijkiidN(0,σ2)where σ2=σb2+σw2K,
K

μiy¯ij

In den Antworten dort (wenn ich richtig verstanden habe) habe ich gelernt, dass die Varianz innerhalb des Subjekts keine Rückschlüsse auf Gruppenmittelwerte hat, und es ist in Ordnung, einfach die Durchschnittswerte der Mittelwerte zur Berechnung des Gruppenmittelwerts zu verwenden, dann die Varianz innerhalb der Gruppe zu berechnen und diese zu verwenden Signifikanztests durchführen. Ich möchte eine Methode verwenden, bei der je größer die Varianz innerhalb des Subjekts ist, desto weniger sicher bin ich, was die Gruppenmittel angeht oder warum es keinen Sinn macht, dies zu wünschen.

σw2

Stéphane Laurent
quelle
Dieses Modell ist sinnvoll, aber der Code scheint diese Informationen nicht zu enthalten, und ich denke, das würde ich gerne verstehen. Mache ich das richtig für die Parameterschätzung?: lmer(Value~Group -1 + (1|Subject), dat) lmer(Value~Group -1 + (1|Subject), dat2)Wobei dat die Originaldaten sind und dat2 die mit kleiner Varianz innerhalb des Subjekts simulierte. Ich bekomme die gleichen Standardfehler.
Flasche
Ich habe es nicht versucht, aber das klingt seltsam. Sie entfernen den festen Abschnitt, aber es gibt einen zufälligen Abschnitt nach Thema. Aus theoretischer Sicht sehe ich kein Problem, aber ich weiß nicht genau, wie ich lmermit Modellen ohne Unterbrechung umgehen soll. Behalten Sie den Abschnitt, um sicherzugehen.
Stéphane Laurent
Ich habe diese Anweisung befolgt , da ich sonst nicht herausfinden konnte, wie ich eine Intervallschätzung erhalten kann. Mein Verständnis der Syntax der R-Formel ist gering, daher macht es möglicherweise keinen Sinn.
Flasche
@Flask AFAIK Derzeit gibt es in R kein Paket, mit dem "korrekte" Konfidenzintervalle für lmerModelle ermittelt werden können. Für Ihr Modell gibt es im speziellen Fall eines ausgewogenen Designs einige exakte Methoden der kleinsten Quadrate, aber ich weiß nicht, ob sie in einem Paket verfügbar sind.
Stéphane Laurent
1
Ich frage mich jedoch, ob das lsmeansPaket zusammen mit dem pbkrtestPaket gute Konfidenzintervalle bieten könnte.
Stéphane Laurent