[R-sig-ME] difficulty fitting a model with nlme()

[Anne Hoen]

Hi, 

I'm working on fitting a model with nlme. I have done this before (with some very helpful advice from this forum and from the excellent Pinheiro & Bates book) on a similar dataset but I'm having troubles with a new set of data. 

I have outcomes measured on test and reference materials made over a series of dilutions and performed over a number of replicates. I want to fit a 3-parameter logistic function to the data (outcome vs dilution) and estimate the effect of test vs reference on the parameters with the replicate number included as a random effect.   

Here is my reproducible example:

library(nlme)

#Data:

x<-c(1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8)
y<-c(1.9302400,1.4943900,1.1843600,1.0083900,0.7697200,0.4876400,0.3005200,0.1877900,1.8066200,1.6685600,1.2516900,1.0409400,0.9356900,0.5950900,0.3552200,0.2214900,1.8101575,1.5273775,1.3062875,1.0022375,0.7663775,0.5120175,0.2972375,0.1930575,1.8740375,1.6142875,1.4112575,1.1021575,0.9438775,0.6029375,0.3582775,0.2234575,1.8090425,1.5115425,1.3347625,0.8812225,0.7459425,0.5252425,0.3252925,0.1809025,1.8420925,1.6119025,1.4094025,1.0518425,0.7759025,0.5792625,0.3502225,0.2061425,1.7288700,1.3824700,1.1691500,0.9555800,0.7044500,0.4966700,0.2938700,0.1374800,1.7447500,1.5038300,1.2124100,0.8477100,0.7546500,0.5371800,0.3137300,0.1925500,2.0704300,1.9151700,1.5235800,1.2380100,0.9800500,0.6039100,0.3916100,0.2459700,2.1507500,1.9126500,1.6317700,1.3445500,0.8570100,0.7079300,0.4359800,0.2923700,2.0717450,1.9443050,1.5423450,1.2771950,0.8258650,0.5431250,0.2090950,0.0916050,2.1291650,1.8801650,1.6445650,1.3246050,0.8991050,0.6972950,0.2931050,0.1032650,1.8966625,1.7421425,1.4608625,1.1705925,0.9664625,0.5481225,0.1529625,0.0761425,2.0282025,1.7920225,1.5809425,1.2019925,0.8361225,0.6382225,0.2437925,0.0803025,2.0854675,1.8130475,1.5290675,1.1757475,0.9437975,0.5506275,0.1783475,0.0795975,2.1994675,1.8947675,1.6205975,1.2414475,0.9887675,0.6347475,0.2842675,0.0842675)
testref<-c("ref","ref","ref","ref","ref","ref","ref","ref","test","test","test","test","test","test","test","test","ref","ref","ref","ref","ref","ref","ref","ref","test","test","test","test","test","test","test","test","ref","ref","ref","ref","ref","ref","ref","ref","test","test","test","test","test","test","test","test","ref","ref","ref","ref","ref","ref","ref","ref","test","test","test","test","test","test","test","test","ref","ref","ref","ref","ref","ref","ref","ref","test","test","test","test","test","test","test","test","ref","ref","ref","ref","ref","ref","ref","ref","test","test","test","test","test","test","test","test","ref","ref","ref","ref","ref","ref","ref","ref","test","test","test","test","test","test","test","test","ref","ref","ref","ref","ref","ref","ref","ref","test","test","test","test","test","test","test","test")
rep<-c("R1","R1","R1","R1","R1","R1","R1","R1","R1","R1","R1","R1","R1","R1","R1","R1","R2","R2","R2","R2","R2","R2","R2","R2","R2","R2","R2","R2","R2","R2","R2","R2","R3","R3","R3","R3","R3","R3","R3","R3","R3","R3","R3","R3","R3","R3","R3","R3","R4","R4","R4","R4","R4","R4","R4","R4","R4","R4","R4","R4","R4","R4","R4","R4","R5","R5","R5","R5","R5","R5","R5","R5","R5","R5","R5","R5","R5","R5","R5","R5","R6","R6","R6","R6","R6","R6","R6","R6","R6","R6","R6","R6","R6","R6","R6","R6","R7","R7","R7","R7","R7","R7","R7","R7","R7","R7","R7","R7","R7","R7","R7","R7","R8","R8","R8","R8","R8","R8","R8","R8","R8","R8","R8","R8","R8","R8","R8","R8")

data<-data.frame(x, y, testref, rep)
data<-groupedData(y~x | rep, data=data) #create a groupedData object

#And the model:

init<-getInitial(y~SSlogis(x, Asym, xmid, scal), data) #find some starting values
fit<-nlme(y~SSlogis(x, Asym, xmid, scal), 
          fixed=Asym + xmid + scal ~ testref-1, #reference vs test as a fixed effect
          random=list(rep=pdDiag(Asym+xmid+scal~1)), #random effect of replicate 
          data=data,
          start=c(init[1], 1, init[2], 1, init[3],1)
)


Interestingly, I don't get an error when I specify the fixed effect as: fixed=list(Asym ~ testref - 1, xmid ~ testref, scal ~ testref); however, this doesn't allow me to estimate the xmid and scal parameters for both the test and reference. i.e.:  
  
fit<-nlme(y~SSlogis(x, Asym, xmid, scal), 
            fixed=list(Asym~testref-1, xmid~testref, scal~testref), #reference vs test as a fixed effect
            random=list(rep=pdDiag(Asym+xmid+scal~1)), #random effect of replicate 
            data=data,
            start=c(init[1], 1, init[2], 1, init[3],1)
  )
summary(fit)
coef(fit)

Nor do I get an error when I specify the original model but remove all observations where x=1, i.e.:
  
data<-data[data$x!=1,]
init<-getInitial(y~SSlogis(x, Asym, xmid, scal), data) #find some starting values
fit<-nlme(y~SSlogis(x, Asym, xmid, scal), 
          fixed=Asym + xmid + scal ~ testref-1, #reference vs test as a fixed effect
          random=list(rep=pdDiag(Asym+xmid+scal~1)), #random effect of replicate 
          data=data,
          start=c(init[1], 1, init[2], 1, init[3],1)
)

summary(fit)
coef(fit)


I am confused. Is this a convergence problem? Is my data not suited to this modeling approach?

I am so stuck! Any help would be very appreciated.