Hi Jenny, (I use the data you provide in the previous e-mail)
For the 1st question, let me assume you only want to compare loc: A vs. B
So you could specified your code like this:
fmAB <- nlme(Y ~ SSlogis(X, Asym, R0, lrc),data = LAST,
random = Asym ~1,
fixed = Asym+R0+lrc ~ loc,
start=c(0.97,0,
1.14,0,
-0.18,0))
summary(fmAB)
Nonlinear mixed-effects model fit by maximum likelihood
Model: Y ~ SSlogis(X, Asym, R0, lrc)
Data: LAST
AIC BIC logLik
-31.02303 -19.70017 24.51152
Random effects:
Formula: Asym ~ 1 | loc
Asym.(Intercept)
StdDev: 1.549779e-06
Formula: Asym ~ 1 | dir %in% loc
Asym.(Intercept) Residual
StdDev: 6.124237e-08 0.09426086
Fixed effects: Asym + R0 + lrc ~ loc
Value Std.Error DF t-value p-value
Asym.(Intercept) 0.9477404 0.1037337 14 9.136280 0.0000
Asym.locB 0.0982456 0.4175971 14 0.235264 0.8174
R0.(Intercept) 1.1289387 0.0652189 14 17.310003 0.0000
R0.locB 0.1390656 0.3946717 14 0.352358 0.7298
lrc.(Intercept) -0.2110057 0.0656513 14 -3.214036 0.0062
lrc.locB -0.0820484 0.6826382 14 -0.120193 0.9060
Then you know Asym, R0, and lrc of loc B are not significant.
Moreover, you can test the joint fixed effect by
anova(fmAB)(Pinherio and Bate, 2000 Book, p 374)
for the 2nd question, How to get the fitted value for particular level?
Based on this example, let me assume you want to get the fitted value of A/N.
then you could write a small code like this:
FV<-data.frame(F.V=fitted(fmAB), group=summary(fmAB)$groups$dir)
A.N<-FV[is.element(FV$group, c("A/N")),]
F.V group
9 0.4209011 A/N
10 0.2726129 A/N
hope this is helpful~
Chunhao
Quoting Jenny Sun <[EMAIL PROTECTED]>:
Thank you for your reply Chunhao!
I attached only part of the test data and that is why you might not
be able to get convergence. Sorry.
I have a couple more questions:
For the second question you answered, how to specify the correct
length of starting values. I tried using the length of levels in
each of the parameters in the start list but found:
fm1 <- nlme(DIFN ~ SSlogis(SVA, Asym, R0, lrc),data = LAST,fixed =
Asym + R0 + lrc ~ dir %in% loc,random = Asym ~ 1,start =list(Asym =
c(1,1,1,1), R0 = c(1,1,1,1), lrc = c(-5,-2,-2,-2)))
Error in nlme.formula(DIFN ~ SSlogis(SVA, Asym, R0, lrc), data = LAST, :
start must have a component called "fixed"
I've got two loc levels (A,B) with four group levels(N,E,S,W); How
I am gonna define the list and the component called"fixed"?
My another question is about the fitted value of the model. If I
want to calculate adjusted R square, I have to get fitted(fm1).
WHich has values like this;
>fitted(fm1)[1:40]
AB/N AB/N AB/E AB/S AB/W AB/W
AB/W AB/W AB/W AB/W
0.6541876 0.7421748 0.8408251 0.5879220 0.4889387 0.6129576
0.5097593 0.6195679 0.5152567 0.5680860
AB/W AB/W AB/W AB/W AB/W AB/W
AB/N AB/N AB/N AB/E
0.4724423 0.8128148 0.7674529 0.7106698 0.6553155 0.6074771
0.5036201 0.5464105 0.6062978 0.6878438
AB/N AB/N AB/N AB/S AB/S AB/S
AB/S AB/S AB/S AB/S
0.7792725 0.8411961 0.7942503 0.7354845 0.5895700 0.6781973
0.6286886 0.5212052 0.8864748 0.8370021
AB/S AB/S AB/N AB/N AB/N AB/N
AB/E AB/E AB/E AB/E
0.7750731 0.7147024 0.6625288 0.5492599 0.5959280 0.6612426
0.7501786 0.8498928 0.6274681 0.7118615
My question is how to get the fitted values for specified group
levels (eg. values for AB/E)?
Thank all very much!
Jenny
Hi Jenny,
I try your code but I did not get in converge in fm3 (see the below).
For the first question, you could use fm1 to interpret the result
without bothering fm2 and fm3. It means that R0 and lrc can be treated
as pure fixed effects (Pinherir and Bates, 2000 Book).
For the second question, your want to know "is AB/E different from the AB/S"
The simplest way is to change your fixed statement:
fixed = Asym+R0+lrc ~ dir %in% loc
and specify the correct length of starting values.
If I am wrong please correct me~
Hope this helpful.
Chunhao Tu
test<-read.table(file="C:\\Documents and
Settings\\ado_cabgfaculty\\Desktop\\sun.txt", header=T)
LAST<-groupedData(Y~X|loc/dir, data=test)
fm1 <- nlme(Y ~ SSlogis(X, Asym, R0, lrc),data = LAST,
+ random = Asym ~1,
+ fixed = Asym+R0+lrc ~ 1,
+ start=c(Asym = 0.97, R0 = 1.14, lrc = -0.18))
fm2 <- update(fm1, random = pdDiag(Asym + R0 ~ 1))
fm3 <- update(fm2, random = pdDiag(Asym+R0+lrc~ 1))
Error in nlme.formula(model = Y ~ SSlogis(X, Asym, R0, lrc), data = LAST, :
Step halving factor reduced below minimum in PNLS step
Quoting Jenny Sun <[EMAIL PROTECTED]>:
My special thanks to Chunhao Tu for the suggestions about testing
significance of two locations.
I used logistic models to describe relationships between Y and X at
two locations (A & B). And within each location, I have four groups
(N,E,S,W)representing directions. So the test data can be arranged as:
Y X dir loc
0.6295 0.8667596 S A
0.7890 0.7324820 S A
0.4735 0.9688875 S A
0.7805 1.1125239 S A
0.8640 0.9506174 E A
0.9445 0.6582157 E A
0.8455 0.5558860 E A
0.9380 0.3304870 E A
0.4010 1.1763090 N A
0.2585 1.3202890 N A
0.3750 1.1763090 E A
0.3855 1.3202890 E A
0.3020 1.1763090 S A
0.2300 1.3202890 S A
0.3155 1.1763090 W A
0.8890 0.6915861 W B
0.9185 0.6149019 W B
0.9275 0.5289258 W B
0.8365 0.9507088 S B
0.7720 0.8842165 N B
0.8615 0.8245123 N B
0.9170 0.7559687 W B
0.9590 0.6772720 W B
0.9900 0.5872023 W B
0.9940 0.4849064 W B
0.7500 0.9560776 W B
The data is grouped using:
LAST<-groupedData(Y~X|loc/dir, data=test)
I then used logistic models to define the relationship between Y and
X, and got fm1, fm2, and fm3 as follows:
--------------------------
fm1 <- nlme(DIFN ~ SSlogis(SVA, Asym, R0, lrc),data = LAST,fixed =
Asym + R0 + lrc ~ 1,random = Asym ~ 1,start =c(Asym = 1, R0 = 1,
lrc = -5))
fm2 <- update(fm1, random = pdDiag(Asym + R0 ~ 1))
fm3 <- update(fm2, random = pdDiag(Asym + R0 + lrc ~ 1))
anova(fm1,fm2,fm3)
------------------------------------------------------------
ANOVA showed:
anova(fm1,fm2,fm3)
Model df AIC BIC logLik Test L.Ratio p-value
fm1 1 7 -1809.913 -1774.304 910.9564
fm2 2 9 -1805.774 -1758.295 910.8871 1 vs 2 0.1386696 0.9999
fm3 3 12 -1801.822 -1742.473 910.9109 2 vs 3 0.0475543 0.9666
** question: do the results show that fm1 could represent the
results of fm2 and fm3?
coef(fm1)
Asym R0 lrc
AB/E 0.9148927 1.389432 -0.3009858
AB/N 0.8775250 1.389432 -0.3009858
AB/S 0.9247592 1.389432 -0.3009858
AB/W 0.8479180 1.389432 -0.3009858
BC/E 0.8791908 1.389432 -0.3009858
BC/N 0.8414229 1.389432 -0.3009858
BC/S 0.9169323 1.389432 -0.3009858
BC/W 0.8817838 1.389432 -0.3009858
** question: how could I know if any of the models is significantly
different from the other ones? (eg. AB/E is different from the AB/S)?
summary(fm1)
Nonlinear mixed-effects model fit by maximum likelihood
Model: DIFN ~ SSlogis(SVA, Asym, R0, lrc)
Data: LAST
AIC BIC logLik
-1809.913 -1774.304 910.9564
Random effects:
Formula: Asym ~ 1 | loc
Asym
StdDev: 2.303402e-05
Formula: Asym ~ 1 | dir %in% loc
Asym Residual
StdDev: 0.03208693 0.1741559
Fixed effects: Asym + R0 + lrc ~ 1
Value Std.Error DF t-value p-value
Asym 0.8855531 0.015375906 2783 57.59355 0
R0 1.3894322 0.009418047 2783 147.52869 0
lrc -0.3009858 0.012833066 2783 -23.45393 0
Correlation:
Asym R0
R0 -0.440
lrc -0.452 0.150
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.1326757 -0.6117037 0.1082112 0.6575250 3.3297270
Number of Observations: 2793
Number of Groups:
loc dir %in% loc
2 8
I have marked all the codes and questions(**). Any answers and
suggestions are appreciated.
Have a good day!
Jenny
______________________________________________
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.
