Dear all,

I have been asked to run an analysis on a dataset with a high number
variables, both continuous and categorical. The data structure, for which I
have included an example matrix below, is balanced, and consists in counts
of bacterial infections (and other non-count variables not included here)
and count of treatment occurrence (DliCount), defined by patient case-match
pair ID (PairID), group of treatment (Dli), and time of measurment (time).
Please find my code and associated output enclosed below.

I am trying to determine which variables (for example, bacterial inections)
significantly differ between groups of treatment (Dli 1 or 0), accounting
for case-matching and repeated measures.

Would you please be kind enough to enlight me about those questions:

Regarding the following equation:
lmer(DliCount~time+Dli+Dli*time+(1|PairID)+(0+inf_bacteria|PairID),data=DLIf,
REML=TRUE, family=gaussian)


   1. Should I state my random effect for subjects as nested (1|PairID/Dli)
   or simply as (1|PairID)? (i.e. Do I correctly assume that nesting is
   already implied?)
   2. When adding random effects that differ among subjects but not in
   intercept (which I assume is 0 for all counts), is the form
   (0+inf_bacteria|PairID) correct?
   3. Does the above equation effectively account for both case-matching
   (on the base of "Dli") and repeated measures ("time")? I am having trouble
   with that since case-matching could be interpreted as a second hierarchical
   level of repeated measures...


Regarding MCMCsamp, I found on google that it cannot define "confidence
intervals" for correlated random effects, but here my random effects are
uncorrelated and mcmcsamp outputs no error but no confidence interval for
inf_bacteria either. How can I obtain those CI?

Thank you very much for reading,
Raoul

CODE+++++++++++++++++++++++++++++++++++++++++++++++++++++++
#packages
require(lme4)

#Example data structure
exampleMatrix<-matrix(c(1,1,1,101,201, 1,1,2,102,202, 1,1,3,103,203,
1,0,1,104,204, 1,0,2,105,205, 1,0,3,106,206,
                2,1,1,107,207, 2,1,2,108,208, 2,1,3,109,209, 2,0,1,110,210,
2,0,2,111,211, 2,0,3,112,212),nrow=12,ncol=5,
    byrow=TRUE,

dimnames=list(c("row1","row2","row3","row4","row5","row6","row7","row8","row9","row10","row11","row12"),
c("PairID","Dli","time","DliCount","inf_bacteria")))
exampleMatrix

#LEGEND
#PairID        ->ID of the patient pair
#Dli            ->Treatment group inside matched pairs
#time            ->time of measurement
#DliCount        ->count of treatment occurrences
#inf_bacteria    ->count of bacterial infections having occurred
#******************************************************************
rm(list=ls()) #removes all objects

#Dataframe constitution
DLI<-read.table("C:\\Users\\Raoul\\Desktop\\DLI
definitif\\donnees_brutes\\multilevel\\outcomeDliCount_2years.csv",
header=TRUE, sep=",")
DLIdf<-data.frame(DLI)
DLIf<-within(DLIdf,{
    PairID<-factor(PairID)
    Dli<-factor(Dli)
    time<-factor(time)
    patID<-interaction(PairID,Dli)
})
str(DLIf)

#Linear mixed model analysis with covariates on DliCount outcome
lmm<-lmer(DliCount~time+Dli+Dli*time+(1|PairID)+(0+inf_bacteria|PairID),data=DLIf,
REML=TRUE, family=gaussian)
print(lmm)
residuals<-resid(lmm)
summary(residuals)
hist(residuals)

#Intraclass correlation for lmm
lmm.null<-lmer(DliCount~1+(1|PairID),data=DLIf)
print(lmm.null)

lmm.null2<-lmer(DliCount~1+(0+inf_bacteria|PairID),data=DLIf)
print(lmm.null2)

#95 CI by MCMC
mcmc.5000<-mcmcsamp(lmm, n=5000)
str(mcmc.5000)
deviance<-as.vector(mcmc.5000@deviance)
summary(deviance) #the smaller the better
HPDinterval(mcmc.5000, prob=0.95)

OUTPUT+++++++++++++++++++++++++++++++++++++++++++++++++++++
 #packages
> require(lme4)
>
> #Example data structure
> exampleMatrix<-matrix(c(1,1,1,101,201, 1,1,2,102,202, 1,1,3,103,203,
1,0,1,104,204, 1,0,2,105,205, 1,0,3,106,206,
+ 2,1,1,107,207, 2,1,2,108,208, 2,1,3,109,209, 2,0,1,110,210,
2,0,2,111,211, 2,0,3,112,212),nrow=12,ncol=5,
+ byrow=TRUE,
+
dimnames=list(c("row1","row2","row3","row4","row5","row6","row7","row8","row9","row10","row11","row12"),
c("PairID","Dli","time","DliCount","inf_bacteria")))
> exampleMatrix
      PairID Dli time DliCount inf_bacteria
row1       1   1    1      101          201
row2       1   1    2      102          202
row3       1   1    3      103          203
row4       1   0    1      104          204
row5       1   0    2      105          205
row6       1   0    3      106          206
row7       2   1    1      107          207
row8       2   1    2      108          208
row9       2   1    3      109          209
row10      2   0    1      110          210
row11      2   0    2      111          211
row12      2   0    3      112          212
>
> #LEGEND
> #PairID->ID of the patient pair
> #Dli->Treatment group inside matched pairs
> #time->time of measurement
> #DliCount->count of treatment occurrences
> #inf_bacteria->count of bacterial infections having occurred
> #******************************************************************
> rm(list=ls()) #removes all objects
>
> #Dataframe constitution
> DLI<-read.table("C:\\Users\\Raoul\\Desktop\\DLI
definitif\\donnees_brutes\\multilevel\\outcomeDliCount_2years.csv",
header=TRUE, sep=",")
> DLIdf<-data.frame(DLI)
> DLIf<-within(DLIdf,{
+ PairID<-factor(PairID)
+ Dli<-factor(Dli)
+ time<-factor(time)
+ patID<-interaction(PairID,Dli)
+ })
> str(DLIf)
'data.frame':   336 obs. of  27 variables:
 $ PairID           : Factor w/ 56 levels "1","2","3","4",..: 1 1 1 2 2 2 3
3 3 4 ...
 $ time             : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1
...
 $ Dli              : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
 $ DliCount         : int  2 2 2 0 2 3 0 0 2 0 ...
 $ TCellDose        : int  100000000 100000000 100000000 0 6000000 16000000
0 0 11000000 0 ...
 $ chimGcas         : num  1 1 1 5.71 4.27 ...
 $ chimMcas         : num  30.2 30.2 30.2 16.9 17.4 ...
 $ CmvDnaCas        : num  308.1 307 NA 57.5 52.5 ...
 $ EbvDna           : num  194 194 NA 0 0 ...
 $ Hemoglobin_mean  : num  96.8 123.5 129.2 97.5 141 ...
 $ WBC_mean         : num  2.93 5.2 5.32 1.65 3.9 ...
 $ NS               : num  0.1111 0.1575 0.1777 0.0165 0 ...
 $ S                : num  1.924 2.921 2.842 0.759 1.716 ...
 $ CD3              : num  0.422 1.437 1.629 0.363 1.792 ...
 $ Monocytes        : num  0.262 0.403 0.436 0.208 0.386 ...
 $ Thrombocytes_mean: num  92.5 175 188.5 141 181 ...
 $ CD4_mean         : num  123 121 315 163 244 ...
 $ CD8_mean         : num  371 352 543 511 406 ...
 $ IgG_mean         : num  10.01 9.29 8.97 7.64 11.8 ...
 $ IgA_mean         : num  1.054 0.949 0.964 1.09 1.3 ...
 $ IgM_mean         : num  0.815 0.853 0.837 1.18 1.05 ...
 $ infections       : int  11 22 NA 1 2 3 3 6 10 2 ...
 $ virus            : int  3 6 NA 1 2 3 0 0 0 0 ...
 $ fungi            : int  3 6 NA 0 0 0 1 2 3 1 ...
 $ inf_bacteria     : int  5 10 NA 0 0 0 2 4 6 1 ...
 $ Neutrophils_total: num  2.036 3.078 3.02 0.775 1.716 ...
 $ patID            : Factor w/ 112 levels "1.0","2.0","3.0",..: 57 57 57
58 58 58 59 59 59 60 ...
>
> #Linear mixed model analysis with covariates on DliCount outcome
>
lmm<-lmer(DliCount~time+Dli+Dli*time+(1|PairID)+(0+inf_bacteria|PairID),data=DLIf,
REML=TRUE, family=gaussian)
> print(lmm)
Linear mixed model fit by REML
Formula: DliCount ~ time + Dli + Dli * time + (1 | PairID) + (0 +
inf_bacteria |      PairID)
   Data: DLIf
   AIC   BIC logLik deviance REMLdev
 737.1 771.3 -359.5    701.6   719.1
Random effects:
 Groups   Name         Variance Std.Dev.
 PairID   (Intercept)  0.12796  0.357713
 PairID   inf_bacteria 0.00170  0.041231
 Residual              0.38241  0.618395
Number of obs: 332, groups: PairID, 56

Fixed effects:
              Estimate Std. Error t value
(Intercept) -0.0008587  0.0958371  -0.009
time2        0.0051175  0.1173416   0.044
time3        0.0245244  0.1200014   0.204
Dli1         0.6869632  0.1173712   5.853
time2:Dli1   0.5874607  0.1657914   3.543
time3:Dli1   1.0977871  0.1691283   6.491

Correlation of Fixed Effects:
           (Intr) time2  time3  Dli1   tm2:D1
time2      -0.600
time3      -0.579  0.501
Dli1       -0.609  0.494  0.480
time2:Dli1  0.427 -0.705 -0.349 -0.695
time3:Dli1  0.416 -0.350 -0.698 -0.674  0.500
> residuals<-resid(lmm)
> summary(residuals)
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-1.6420000 -0.2632000 -0.0341500 -0.0000031  0.2928000  2.5470000
> hist(residuals)
>
> #Intraclass correlation for lmm
> lmm.null<-lmer(DliCount~1+(1|PairID),data=DLIf)
> print(lmm.null)
Linear mixed model fit by REML
Formula: DliCount ~ 1 + (1 | PairID)
   Data: DLIf
   AIC  BIC logLik deviance REMLdev
 996.7 1008 -495.3    986.8   990.7
Random effects:
 Groups   Name        Variance   Std.Dev.
 PairID   (Intercept) 2.9744e-12 1.7246e-06
 Residual             1.1075e+00 1.0524e+00
Number of obs: 336, groups: PairID, 56

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.66369    0.05741   11.56
>
> lmm.null2<-lmer(DliCount~1+(0+inf_bacteria|PairID),data=DLIf)
> print(lmm.null2)
Linear mixed model fit by REML
Formula: DliCount ~ 1 + (0 + inf_bacteria | PairID)
   Data: DLIf
   AIC   BIC logLik deviance REMLdev
 947.5 958.9 -470.8    937.7   941.5
Random effects:
 Groups   Name         Variance Std.Dev.
 PairID   inf_bacteria 0.010705 0.10347
 Residual              0.841557 0.91736
Number of obs: 332, groups: PairID, 56

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.50678    0.05821   8.706
>
> #95 CI by MCMC
> mcmc.5000<-mcmcsamp(lmm, n=5000)
> str(mcmc.5000)
Formal class 'merMCMC' [package "lme4"] with 9 slots
  ..@ Gp      : int [1:3] 0 56 112
  ..@ ST      : num [1:2, 1:5000] 0.5785 0.0667 0.6265 0.0724 0.4446 ...
  ..@ call    : language lmer(formula = DliCount ~ time + Dli + Dli * time
+ (1 | PairID) + (0 +      inf_bacteria | PairID), data = DLIf, REML = TRUE)
  ..@ deviance: num [1:5000] 702 702 702 704 706 ...
  ..@ dims    : Named int [1:18] 2 332 6 112 1 2 1 1 2 5 ...
  .. ..- attr(*, "names")= chr [1:18] "nt" "n" "p" "q" ...
  ..@ fixef   : num [1:6, 1:5000] -0.000859 0.005118 0.024524 0.686963
0.587461 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:6] "(Intercept)" "time2" "time3" "Dli1" ...
  .. .. ..$ : NULL
  ..@ nc      : int [1:2] 1 1
  ..@ ranef   : num[1:112, 0 ]
  ..@ sigma   : num [1, 1:5000] 0.618 0.512 0.591 0.607 0.63 ...
> deviance<-as.vector(mcmc.5000@deviance)
> summary(deviance) #the smaller the better
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  701.6   705.1   707.6   708.3   710.6   730.7
> HPDinterval(mcmc.5000, prob=0.95)
$fixef
                 lower     upper
(Intercept) -0.1805506 0.1820464
time2       -0.2477150 0.2287558
time3       -0.2377523 0.2567687
Dli1         0.4479843 0.9217012
time2:Dli1   0.2516018 0.9374890
time3:Dli1   0.7637042 1.4634564
attr(,"Probability")
[1] 0.95

$ST
          lower     upper
[1,] 0.21267994 0.5363795
[2,] 0.01438342 0.1143554
attr(,"Probability")
[1] 0.95

$sigma
         lower     upper
[1,] 0.5917736 0.7119486
attr(,"Probability")
[1] 0.95

	[[alternative HTML version deleted]]

