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:
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:
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:
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:
Many thanks in advance for any help.
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)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)
Code:
> summary(glht(excit.model.1, linfct = mcp(treat.f = "Tukey")))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 ctrlMany thanks in advance for any help.