Gastric Emptying Analysis with R: Miscellaneous

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 write

nlsList(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 print at the package author's discretion. When we plan to do some more work with the data, this is the way to go:

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.998
Residual  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² to show the quality of the fit!

To quote Douglas Bates:

... "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.