0
$\begingroup$

We have done an experiment using optogenetics (a technique to manipulate genetically engineered neurons with light) and are trying to write the proper mixed-effects model.

The experiment is as follows:

  • subjects are assigned to one of two tasks (A or E)
  • on each day the subject undergoes a single session where the the fiber-optic cable is connected region X or Y
  • In a single session, there are trials with laser on and off (1 or 0)

We are interested whether the difference in effect between X and Y are different between the two tasks A and E and we want to account for the fact that performance can vary from session to session, so we want to use the no laser trials to adjust for across session variation. There is also across subject variation (some have overall higher performance than others).

We have tried models like p ~ 0 + task&laser&region + (1|subjid) + (1|sessid) which provides an estimate of the task/laser/region effect accounting for variation across sessions and subject. Then in order to get the statistic we want we have to compare the coefficients for the different interaction terms.

Is this the right model? Seems it is missing something.

Example csv and the code below to generate it.

subjid,sessid,task,region,laser,p       ,n
     2,    93,A   ,Y     ,    0,0.444444,18
     2,    93,A   ,Y     ,    1,0.285714,21
     2,    99,A   ,X     ,    0,0.684211,19
     2,    99,A   ,X     ,    1,0.607143,28
     4,    91,A   ,Y     ,    0,0.8     ,30
     4,    91,A   ,Y     ,    1,0.538462,26
     4,    97,A   ,X     ,    0,0.85    ,20
     4,    97,A   ,X     ,    1,0.761905,21
     6,    89,A   ,Y     ,    0,0.888889,18
     6,    89,A   ,Y     ,    1,0.7     ,30
     6,    95,A   ,X     ,    0,1.0     ,12
     6,    95,A   ,X     ,    1,0.84    ,25
     1,    94,E   ,Y     ,    0,0.72    ,25
     1,    94,E   ,Y     ,    1,0.307692,13
     1,   100,E   ,X     ,    0,0.6     ,10
     1,   100,E   ,X     ,    1,0.5     ,10
     3,    92,E   ,Y     ,    0,0.727273,11
     3,    92,E   ,Y     ,    1,0.434783,23
     3,    98,E   ,X     ,    0,1.0     ,12
     3,    98,E   ,X     ,    1,0.9375  ,16
     5,    90,E   ,Y     ,    0,0.875   ,16
     5,    90,E   ,Y     ,    1,0.933333,15
     5,    96,E   ,X     ,    0,0.956522,23
     5,    96,E   ,X     ,    1,0.965517,29
jdf = let
   task = ["A", "E"]
   laser = [0,1]
   fiber = ["X", "Y"]
   subjid = 1:6
   sessids =    collect(1:100)
    D = Dict{Symbol,Any}[]
    _f(x) = begin
        subjid, fiber = x
        sessid = pop!(sessids) 
        [Dict(
            :subjid => subjid,
            :task => task[(subjid % 2 + 1)],
            :laser => l,
            :region => fiber,
            :sessid => sessid,
            :n => rand(10:30),
            :gen_p => logit(randn()*2 + 2)
            ) for l in laser]
    end
    todo = [(i,j) for i in subjid, j in fiber]
    foreach((x)->push!(D,_f(x)...), todo)
    @chain DataFrame(D) begin
        @rtransform(:p = rand(Binomial(:n,:gen_p))/:n)
        select(
        :subjid, :sessid, :task, :region, :laser, :p, :n, :gen_p)
    end
end
$\endgroup$
1
  • $\begingroup$ Reproducibility is great but in this case it may be more productive to paste the generated data instead the julia code. $\endgroup$
    – dipetkov
    Commented Apr 30 at 19:54

1 Answer 1

0
$\begingroup$

To test different models (and i tested many!) i generated synthetic data with known interactions and then compared the AIC across models.

This helped me realize that the correct thing was to compare the following nested models:

p ~ 0 + task + laser  + laser&region + laser:region:task + (1|subjid/sessid)
p ~ 0 + task + laser  + laser&region + (1|subjid/sessid)

Both models have terms to capture overall task differences, the effect of laser and an interaction between laser and region. Only the first has an extra interaction term for the 3-way interaction of laser:region:task.

Then a likelihood ratio test

Likelihood-ratio test: 
───────────────────────────────────────────────────────
     DOF  ΔDOF     LogLik  Deviance    Chisq  p(>Chisq)
───────────────────────────────────────────────────────
[1]    6        -110.5747   81.7305                    
[2]    8     2   -93.9950   48.5711  33.1593     <1e-07
───────────────────────────────────────────────────────
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.