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
______________________________________________
R-help@r-project.org 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.
