### read.table separators

When your data columns are separated by tabs instead of commas, use:

ge = read.table("gastempt.csv",sep="\t",header=TRUE)

### Reading from a database

With package
RODBC,
you can read in data from Excel, Access or any other database for which an ODBC driver is installed. For example, to read in data from a range named *data* in an Excel file, use:

library(RODBC) channel = odbcConnectExcel("haircount.xls") hc = sqlQuery(channel,"select * from data") odbcClose(channel)

We recommend using named ranges for reading from Excel, reading worksheets in theory is possible, but does not always work as expected.

### When to assign and when to print.

When you want to compte a fit with nlsList, you could simply writenlsList(v~SSEmptLinExp(t,v0,kappa,tempt)|subjtreat,data=ge)

When the result is not assigned to a left-side value by "<-" or "=", it is implicitely assumed that the function is handled by some internal print function which knows how to handle the type of output. So will will get the estimate table printed, and that´s it, we cannot do more work with the result. The above thus is equivalent to

print(nlsList(v~SSEmptLinExp(t,v0,kappa,tempt)|subjtreat,data=ge))If you want some more output, you can write

summary(nlsList(v~SSEmptLinExp(t,v0,kappa,tempt)|subjtreat,data=ge))The set of

*summary*functions always give a bit more than

gelist <- nlsList(v~SSEmptLinExp(t,v0,kappa,tempt)|subjtreat,data=ge) gelist = nlsList(v~SSEmptLinExp(t,v0,kappa,tempt)|subjtreat,data=ge)The first form is prefered by dyed-in-the wool R oldtimers; I prefer the second form, since it looks more like other programming languages. We can use the result assigned to

*gelist*now for further processing, but nothing is printed; when we assign, we have to explicitely request a printout by one of the following methods

gelist print(gelist) summary(gelist)

The first form, simply quoting the variable, is most common in interactive mode. The second gives the same output, but make more explicit what happens and is sometimes required. The third form gives extended information.

So here comes the form we used in the program:

gelist = nlsList(v~SSEmptLinExp(t,v0,kappa,tempt)|subjtreat,data=ge) print(gelist)

### Controlling nlme convergence

In a perfect world, computing an *nlme* fit from a given result of *nlsList* should be as easy as

ge0.nlme=nlme(gelist)

The real world requires some tweaking. The default value 0.001 of parameter *pnlsTol*
controling the first steps of the nlme fit was too small for our data set,
and lead to oscillations in convergence
(ping-ponging according to Douglas Bates. Try *verbose=T* in case of non-convergence to see the missed approach, and use smaller values of *pnlsTol* if possible.

contr = nlmeControl(pnlsTol=0.3) ge0.nlme=nlme(gelist,control=contr,verbose=T)

### Correlation of Coefficients

The result of fitting Example 1 with nlme are not fully satisfactory, if we look a bit closer. Here is the summary output of the fit, heavily shortened:

> summary(ge0.nlme) Nonlinear mixed-effects model fit by maximum likelihood Model: v ~ SSEmptLinExp(t, v0, kappa, tempt) Random effects: Formula: list(v0 ~ 1, kappa ~ 1, tempt ~ 1) StdDev Corr v0 189.664 v0 kappa kappa 0.678 -0.247 tempt 39.624 0.181-0.998Residual 29.133 ... Fixed effects: list(v0 ~ 1, kappa ~ 1, tempt ~ 1) Value Std.Error DF t-value p-value v0 644 68.6 78 9.38 0 kappa 1 0.3 78 5.70 0 tempt 75 14.6 78 5.11 0 ... Number of Observations: 88 Number of Groups: 8

The fixed effect coefficients are secondary here, because these are group averages and we are chiefly interested in the coefficients of the single fit. However, the correlation coefficients of the random effect (and of the fixed effects, not shown here) show a very high correlation of 0.998 between kappa and t_{emtp}. Such a correlation indicates that somehow we tried to fit too many things, or that our emptying curves lack personality required to make three parameters happy; See PB chapter 8.2 for a discussion and and treatment.

This is partially also due to the selection we made of the emptying data in the study; for the fit of our full 72 emptying curves, the correlation was much lower, because there were more curves with the initial rise that parameter *kappa* models. The essential idea to improve the result is that we don't have 8 groups in our data set, as the last line suggests, but 2 treatments and 4 subjects. We expect that subjects and treatments do have something in common, and we should tell *nlme* about this. Thus, to get good estimates, we must put the concept of or study into the fit. This is not easy, as you will certainly see when you have read PB chapter 8.2, but we will give it a try in a later example.

### You should use R^{2} to show the quality of the fit!

- ... "the lack of automatic ANOVA, R^2 and adj. R^2 from nls is a feature, not a bug :-)"
- My best advice regarding R^2 statistics with nonlinear models is, as Nancy Reagan suggested, "Just say no."

To avoid further problems with punctuation, the second quotation has not been quoted.