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®ion + (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