[R] Setting random effects within a category using nlme

[Michael A. Gilchrist]

Hello,

I will start out with the caveat that I'm not a statistician by training, but 
have a fairly decent understanding of probability and likelihood.

Nevertheless, I'm trying to fit a nonlinear model to a dataset which has two 
main factors using nlme. Within the dataset there are two Type categories and 
four  Tissue categories, thus giving me 8 datasets in total.  The dataset is 
in a groupedData object with
 	       formula= Count ~ Time|Type/Tissue
and there are three basic model parameters that I am trying to fit: T0, aN, 
and aL.

Calling nlsList gives
------------------------------------------
>nlsList(Count ~ quad.PBMC.model(aL, aN, T0), data = tissueData, start = 
list(T0 = 1000, aN = exp(-2), aL = exp(-2)))
Call:
   Model: Count ~ quad.PBMC.model(aL, aN, T0) | Type/Tissue
    Data: tissueData

Coefficients:
                     T0         aN         aL
Naive/CLN    1530.0088 0.26508876 0.04730525
Naive/ILN     296.4755 0.09270158 0.09542535
Naive/IngLN   828.1406 0.50799864 0.12500593
Naive/MLN     487.6856 0.16565269 0.10385125
Memory/CLN   3567.2132 0.05656948 0.07753467
Memory/ILN    708.1642 0.01264033 0.10018441
Memory/IngLN 2114.1868 0.05298126 0.12795589
Memory/MLN   1050.0811 0.02277018 0.13560749

Degrees of freedom: 96 total; 72 residual
Residual standard error: 271.4194
>
---------------------------------------------

I find that T0 varies greatly with each treatment, so I'm going to leave that 
alone for now. However, aN and aL values don't seem to vary w/in a Type or 
between Types. As a result, I would like to fit two mixed effects models.  The 
first model assumes that aN and aL are from the same population, independent 
of Type.  I believe that this is doing what I want,


------------------------------------------
model1 = nlme(Count ~ quad.PBMC.model(aL, aN, T0),
      data = tissueData,
      start = list( fixed = c(rep(1000, 8), -2, -2) ),
      fixed = list(T0 ~ TypeTissue-1, aL ~ 1, aN ~ 1),
      random = aL + aN ~ 1|Tissue
      )


Nonlinear mixed-effects model fit by maximum likelihood
   Model: Count ~ quad.PBMC.model(aL, aN, T0)
   Data: tissueData
   Log-likelihood: -669.9258
   Fixed: list(T0 ~ TypeTissue - 1, aL ~ 1, aN ~ 1)
   T0.TypeTissueMemory/CLN   T0.TypeTissueMemory/ILN T0.TypeTissueMemory/IngLN
               3517.767851                692.356073               2058.228897
   T0.TypeTissueMemory/MLN    T0.TypeTissueNaive/CLN    T0.TypeTissueNaive/ILN
                998.151121               1677.199160                301.686695
  T0.TypeTissueNaive/IngLN    T0.TypeTissueNaive/MLN                        aL
                841.100309                496.307552                 -2.409191
                        aN
                 -3.695229

Random effects:
  Formula: list(aL ~ 1, aN ~ 1)
  Level: Tissue
  Structure: General positive-definite, Log-Cholesky parametrization
          StdDev       Corr
aL       2.148064e-01 aL
aN       4.548214e-04 -0.001
Residual 2.544745e+02

Number of Observations: 96
Number of Groups: 4

------------------------------------------
Note that I do not know how to get a separate T0 fit for each Type/Tissue, so 
I have created a column TypeTissue that has a unique designator for each 
Type/Tissue combination.  Nevertheless, I *think* this is doing what I want. 
Given my limited stats training, I also like the fact that the Corr b/w aL and 
aN is very low.



Moving on, the second model I would like to fit assumes that the population of 
aN and aL values across the set of Tissues differ between Types.  I have tried 
a number of different syntaxes, but I can't seem to get the output I was 
expecting. For example, if I run

------------------------------------------
> model2 = nlme(Count ~ quad.PBMC.model(aL, aN, T0),
+      data = tissueData,
+      start = list( fixed = c(rep(1000, 8),rep(-2, 2), rep(-2, 2) )),
+      fixed = list(T0 ~ TypeTissue-1, aL ~ Type-1, aN ~ Type-1),
+      random = aL + aN ~ 1|Tissue
+      )
Nonlinear mixed-effects model fit by maximum likelihood
   Model: Count ~ quad.PBMC.model(aL, aN, T0),
   Data: tissueData
   Log-likelihood: -665.1197
   Fixed: list(T0 ~ TypeTissue - 1, aL ~ Type - 1, aN ~ Type - 1)
   T0.TypeTissueMemory/CLN   T0.TypeTissueMemory/ILN T0.TypeTissueMemory/IngLN
               3587.056207                711.544775               2089.597686
   T0.TypeTissueMemory/MLN    T0.TypeTissueNaive/CLN    T0.TypeTissueNaive/ILN
               1027.412736               1554.219211                271.100935
  T0.TypeTissueNaive/IngLN    T0.TypeTissueNaive/MLN              aL.TypeNaive
                782.986653                460.000116                 -2.698370
             aL.TypeMemory              aN.TypeNaive             aN.TypeMemory
                 -2.253364                 -1.718506                 -3.318075

Random effects:
  Formula: list(aL ~ 1, aN ~ 1)
  Level: Tissue
  Structure: General positive-definite, Log-Cholesky parametrization
                StdDev      Corr
aL.(Intercept)   0.1997321 aL.(I)
aN.(Intercept)   0.3077296 -0.994
Residual       240.7163763

Number of Observations: 96
Number of Groups: 4
>
------------------------------------------

I do get separate estimates of aL and aN for each Type. However, I'm at a loss 
as to how to interpret the "Random effects" part which makes me think I'm not 
doing what I wanted.  To be more specific, I was expecting to see an estimate 
for the StdDev of aL.TypeMemory, aLTypeNaive, aN.TypeMemory, and aNTypeNaive, 
but instead I get Intercepts, which confuses me.  Note that the values in the 
Type and Tissue categories in the dataframe are strings, but I haven't 
knowingly designated them as Factors.

If someone could set me straight on this it would be greatly appreciated.

Sincerely,

Mike

-----------------------------------------------------
Department of Ecology & Evolutionary Biology
569 Dabney Hall
University of Tennessee
Knoxville, TN 37996-1610

phone:(865) 974-6453
fax:  (865) 974-6042

web: http://eeb.bio.utk.edu/gilchrist.asp