[R – lme function] Comparing curves in longitudinal data (growth curves)

#1
I think my problem is best solved by the use of multi-level models. However, I’d like to confirm if what I’m doing is right, and the use of the function is right too.
I work with cells (biology), by doing electrophysiology. I apply to each cell a series of currents, in 20 pA steps, ranging from 0 to 200 pA. At each step the cell produces a number of action potential (AP), and the whole sequence is called excitability of the cell. I’d like to see if the application of certain drugs (treatments) has any effect on this excitability (attached is a .jpg image of this).
What I’ve learned from multi-level models is that in this case, I should consider the treatments and the current steps as fixed factors, and the cells tested as random factor.
What I’ve done so far is using the following code:
Code:
> excit.model.1<- lme(APs ~ StepCurrentpA+StepCurrentpA^2+treat.f, random = ~1 | ID.Cell, data = cluster5.excit,method="ML", na.action=na.omit)
Here’s one question: when I try to specify the random model with both intercepts and slope (random =~StepCurrentpA|ID.Cell), the calculation takes forever, and when it finishes, it’s because the iterations ended. Why does that happen? My guess is that somehow having random slope with current steps makes the slopes parallel and it can’t ever compute the intercepts.
I’ve read that one way of comparing curves would be by doing an ANOVA from the previous model to another one. The last model was, as far as I understood, comparing the effects on the intercept. Then I’d have to compare the effects on current steps, and do an ANOVA between the two:
Code:
> excit.model.2<- lme(APs ~ StepCurrentpA+StepCurrentpA^2+treat.f+StepCurrentpA*treat.f, random = ~1 | ID.Cell, data = cluster5.excit,method="ML",na.action=na.omit)
Code:
> anova(excit.model.1,excit.model.2)
This is all good, but I can’t see it showing me that the drugs have an effect on the excitability of the cells. I tried to do a Tukey multiple comparisons of means:
Code:
> summary(glht(excit.model.1, linfct = mcp(treat.f = "Tukey")))
However, the result it gives me is still for an average of the current steps.
Is there a better way to compare this type of curves, or is there anything in my understanding that is totally wrong?
Below are some info about the data:
Code:
> summary(cluster5.excit)
HTML:
ID.Cell	Treatment.Code	StepCurrentpA	APs	treat.f   
760:143	Min.: 1.000	0: 214	Min.: 0.000	ctrl: 881
788:143	1st Qu.: 1.000	20: 214	1st Qu.: 0.000	5HT: 165  
789:143	Median: 4.000	40: 214	Median: 0.000	DEA: 448  
790:143	Mean: 5.224	60: 214	Mean: 2.483	5HT/DEA: 299  
791:143	3rd Qu.: 8.000	80: 214	3rd Qu.: 5.000	ODQ/DEA: 330  
825:132	Max.: 14.000	100: 214	Max.: 15.000	ODQ/5HT/DEA: 231  
(Other):1507	(Other):1070
Code:
 > str(cluster5.excit)
HTML:
data.frame': 2354 obs. of  13 variables:
 $ ID.Cell: Factor w/ 26 levels "434","435","558",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Treatment.Code: int  1 1 1 1 1 1 1 1 1 1 ...
 $ StepCurrentpA: Factor w/ 11 levels "0","20","40",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ APs: int  0 0 0 0 0 0 0 1 5 7 ...
 $ treat.f: Factor w/ 6 levels "ctrl","5HT","DEA",..: 1 1 1 1 1 1 1 1 1 1 ...
Code:
> head(cluster5.excit)
HTML:
ID.Cell	Treatment.Code	StepCurrentpA	APs	treat.f
434	1	0	0	ctrl
434	1	20	0	ctrl
434	1	40	0	ctrl
434	1	60	0	ctrl
434	1	80	0	ctrl
434	1	100	0	ctrl

Many thanks in advance for any help.