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.
