From fa48699d77af7fa8f12366b2ce5ea555bb5fcf19 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 2 Feb 2024 12:03:27 +0000 Subject: [PATCH] JuliaFormatter --- docs/literate/HowTo/effects.jl | 29 +- docs/literate/HowTo/juliacall_unfold.jl | 2 +- docs/literate/HowTo/timesplines.jl | 1 + docs/literate/HowTo/unfold_io.jl | 16 +- .../explanations/nonlinear_effects.jl | 70 +-- docs/literate/explanations/window_length.jl | 129 ++++-- docs/make.jl | 90 ++-- docs/nb_timesplines.jl | 120 +++--- docs/nb_timesplines_dev.jl | 71 +-- docs/nb_unfold-overlap.jl | 70 +-- docs/nb_unfold-overlap_extended.jl | 92 ++-- docs/notebooks/nb_nonlinear.jl | 52 +-- docs/run_liveserver.jl | 13 +- .../UnfoldBSplineKitExt.jl | 21 +- ext/UnfoldBSplineKitExt/basisfunctions.jl | 5 +- ext/UnfoldBSplineKitExt/splinepredictors.jl | 30 +- ext/UnfoldKrylovExt.jl | 126 +++--- .../UnfoldMixedModelsExt.jl | 4 +- ext/UnfoldMixedModelsExt/condense.jl | 16 +- ext/UnfoldMixedModelsExt/designmatrix.jl | 80 ++-- ext/UnfoldMixedModelsExt/effects.jl | 2 +- ext/UnfoldMixedModelsExt/fit.jl | 29 +- ext/UnfoldMixedModelsExt/statistics.jl | 127 +++--- ext/UnfoldMixedModelsExt/timeexpandedterm.jl | 2 +- ext/UnfoldMixedModelsExt/typedefinitions.jl | 14 +- ext/UnfoldRobustModelsExt.jl | 44 +- src/Unfold.jl | 76 ++-- src/basisfunctions.jl | 57 ++- src/condense.jl | 20 +- src/designmatrix.jl | 241 ++++++----- src/effects.jl | 174 ++++---- src/eventhandling.jl | 35 +- src/fit.jl | 80 ++-- src/io.jl | 18 +- src/predict.jl | 110 ++--- src/solver.jl | 36 +- src/splinepredictors.jl | 1 + src/statistics.jl | 1 + src/timeexpandedterm.jl | 10 +- src/typedefinitions.jl | 2 +- src/utilities.jl | 45 +- test/basisfunctions.jl | 44 +- test/designmatrix.jl | 405 +++++++++--------- test/effects.jl | 316 ++++++++------ test/fit.jl | 91 ++-- test/fit_LMM.jl | 128 ++++-- test/io.jl | 105 +++-- test/predict.jl | 16 +- test/runtests.jl | 2 +- test/setup.jl | 4 +- test/splines.jl | 14 +- test/statistics.jl | 69 +-- test/time_expand.jl | 11 +- test/utilities.jl | 40 +- 54 files changed, 1889 insertions(+), 1517 deletions(-) diff --git a/docs/literate/HowTo/effects.jl b/docs/literate/HowTo/effects.jl index 5e8e21b9..1a0c8a06 100644 --- a/docs/literate/HowTo/effects.jl +++ b/docs/literate/HowTo/effects.jl @@ -14,14 +14,14 @@ using UnfoldMakie # # Setup things # Let's generate some data and fit a model of a 2-level categorical and a continuous predictor with interaction. -data,evts = UnfoldSim.predef_eeg(;noiselevel=8) +data, evts = UnfoldSim.predef_eeg(; noiselevel = 8) basisfunction = firbasis(τ = (-0.1, 0.5), sfreq = 100, name = "basisA") -f = @formula 0 ~ 1+condition+continuous # 1 +f = @formula 0 ~ 1 + condition + continuous # 1 -m = fit(UnfoldModel, Dict(Any=>(f,basisfunction)), evts, data,eventcolumn="type") +m = fit(UnfoldModel, Dict(Any => (f, basisfunction)), evts, data, eventcolumn = "type") # Plot the results plot_erp(coeftable(m)) @@ -30,23 +30,28 @@ plot_erp(coeftable(m)) # ### Effects # A convenience function is `effects`. It allows to specify effects on specific levels, while setting non-specified ones to a typical value (usually the mean) -eff = effects(Dict(:condition => ["car","face"]),m) -plot_erp(eff;mapping=(;color=:condition,)) +eff = effects(Dict(:condition => ["car", "face"]), m) +plot_erp(eff; mapping = (; color = :condition,)) # We can also generate continuous predictions -eff = effects(Dict(:continuous => -5:0.5:5),m) -plot_erp(eff;mapping=(;color=:continuous,group=:continuous=>nonnumeric),categorical_color=false,categorical_group=false) +eff = effects(Dict(:continuous => -5:0.5:5), m) +plot_erp( + eff; + mapping = (; color = :continuous, group = :continuous => nonnumeric), + categorical_color = false, + categorical_group = false, +) # or split it up by condition -eff = effects(Dict(:condition=>["car","face"],:continuous => -5:2:5),m) -plot_erp(eff;mapping=(;color=:condition,col=:continuous=>nonnumeric)) +eff = effects(Dict(:condition => ["car", "face"], :continuous => -5:2:5), m) +plot_erp(eff; mapping = (; color = :condition, col = :continuous => nonnumeric)) # ## What is typical anyway? # The user can specify the typical function applied to the covariates/factors that are marginalized over. This offers even greater flexibility. # Note that this is rarely necessary, in most applications the mean will be assumed. -eff_max = effects(Dict(:condition=>["car","face"]),m;typical=maximum) # typical continuous value fixed to 1 +eff_max = effects(Dict(:condition => ["car", "face"]), m; typical = maximum) # typical continuous value fixed to 1 eff_max.typical .= :maximum -eff = effects(Dict(:condition=>["car","face"]),m) +eff = effects(Dict(:condition => ["car", "face"]), m) eff.typical .= :mean # mean is the default -plot_erp(vcat(eff,eff_max);mapping=(;color=:condition,col=:typical)) +plot_erp(vcat(eff, eff_max); mapping = (; color = :condition, col = :typical)) diff --git a/docs/literate/HowTo/juliacall_unfold.jl b/docs/literate/HowTo/juliacall_unfold.jl index ab9e10a9..c27980c8 100644 --- a/docs/literate/HowTo/juliacall_unfold.jl +++ b/docs/literate/HowTo/juliacall_unfold.jl @@ -37,4 +37,4 @@ # ``` # ## Example: Unfold model fitting from Python -# In this [notebook](https://github.com/unfoldtoolbox/Unfold.jl/blob/main/docs/src/HowTo/juliacall_unfold.ipynb), you can find a more detailed example of how to use Unfold from Python to load data, fit an Unfold model and visualise the results in Python. \ No newline at end of file +# In this [notebook](https://github.com/unfoldtoolbox/Unfold.jl/blob/main/docs/src/HowTo/juliacall_unfold.ipynb), you can find a more detailed example of how to use Unfold from Python to load data, fit an Unfold model and visualise the results in Python. diff --git a/docs/literate/HowTo/timesplines.jl b/docs/literate/HowTo/timesplines.jl index e69de29b..8b137891 100644 --- a/docs/literate/HowTo/timesplines.jl +++ b/docs/literate/HowTo/timesplines.jl @@ -0,0 +1 @@ + diff --git a/docs/literate/HowTo/unfold_io.jl b/docs/literate/HowTo/unfold_io.jl index e4fcfb40..f28abec7 100644 --- a/docs/literate/HowTo/unfold_io.jl +++ b/docs/literate/HowTo/unfold_io.jl @@ -10,14 +10,14 @@ # ### Simulate some example data using UnfoldSim.jl using UnfoldSim -data, events = UnfoldSim.predef_eeg(; n_repeats=10) -first(events,5) +data, events = UnfoldSim.predef_eeg(; n_repeats = 10) +first(events, 5) # ### Fit an Unfold model using Unfold -basisfunction = firbasis(τ = (-0.5,1.0), sfreq = 100, name = "stimulus") -f = @formula 0 ~ 1 + condition + continuous -bfDict = Dict(Any => (f,basisfunction)) +basisfunction = firbasis(τ = (-0.5, 1.0), sfreq = 100, name = "stimulus") +f = @formula 0 ~ 1 + condition + continuous +bfDict = Dict(Any => (f, basisfunction)) m = fit(UnfoldModel, bfDict, events, data); # ```@raw html @@ -28,9 +28,9 @@ m = fit(UnfoldModel, bfDict, events, data); # The following code saves the model in a compressed .jld2 file. The default option of the `save` function is `compress=false`. # For memory efficiency the designmatrix is set to missing. If needed, it can be reconstructed when loading the model. -save_path = mktempdir(; cleanup=false) # create a temporary directory for the example -save(joinpath(save_path, "m_compressed.jld2"), m; compress=true); +save_path = mktempdir(; cleanup = false) # create a temporary directory for the example +save(joinpath(save_path, "m_compressed.jld2"), m; compress = true); # The `load` function allows to retrieve the model again. # By default, the designmatrix is reconstructed. If it is not needed set `generate_Xs=false`` which improves time-efficiency. -m_loaded = load(joinpath(save_path, "m_compressed.jld2"), UnfoldModel, generate_Xs=true); \ No newline at end of file +m_loaded = load(joinpath(save_path, "m_compressed.jld2"), UnfoldModel, generate_Xs = true); diff --git a/docs/literate/explanations/nonlinear_effects.jl b/docs/literate/explanations/nonlinear_effects.jl index 1ebfc6fb..61f7c277 100644 --- a/docs/literate/explanations/nonlinear_effects.jl +++ b/docs/literate/explanations/nonlinear_effects.jl @@ -2,7 +2,7 @@ -using BSplineKit,Unfold +using BSplineKit, Unfold using CairoMakie using DataFrames using Random @@ -14,42 +14,42 @@ using Colors # We start with generating rng = MersenneTwister(2) # make repeatable n = 20 # datapoints -evts = DataFrame(:x=>rand(rng,n)) -signal = -(3*(evts.x .-0.5)).^2 .+ 0.5 .* rand(rng,n) +evts = DataFrame(:x => rand(rng, n)) +signal = -(3 * (evts.x .- 0.5)) .^ 2 .+ 0.5 .* rand(rng, n) -plot(evts.x,signal) +plot(evts.x, signal) # # Looks perfectly non-linear - great! # # # Compare linear & non-linear fit # First we have to reshape it to fit to Unfold format, a 3d array, 1 channel x 1 timepoint x 20 datapoints -signal = reshape(signal,length(signal),1,1) -signal = permutedims(signal,[3,2,1]) +signal = reshape(signal, length(signal), 1, 1) +signal = permutedims(signal, [3, 2, 1]) size(signal) # Next we define three different models. **linear** **4 splines** and **10 splines** # note the different formulas, one `x` the other `spl(x,4)` -design_linear = Dict(Any=> (@formula(0~1+x),[0])) -design_spl3 = Dict(Any=> (@formula(0~1+spl(x,4)),[0])) -design_spl10 = Dict(Any=> (@formula(0~1+spl(x,10)),[0])) #hide +design_linear = Dict(Any => (@formula(0 ~ 1 + x), [0])) +design_spl3 = Dict(Any => (@formula(0 ~ 1 + spl(x, 4)), [0])) +design_spl10 = Dict(Any => (@formula(0 ~ 1 + spl(x, 10)), [0])) #hide # and we fit the parameters -uf_linear = fit(UnfoldModel,design_linear,evts,signal); -uf_spl3 = fit(UnfoldModel,design_spl3,evts,signal); -uf_spl10 = fit(UnfoldModel,design_spl10,evts,signal); #hide +uf_linear = fit(UnfoldModel, design_linear, evts, signal); +uf_spl3 = fit(UnfoldModel, design_spl3, evts, signal); +uf_spl10 = fit(UnfoldModel, design_spl10, evts, signal); #hide # extract the fitted values using Unfold.effects -p_linear = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_linear); -p_spl3 = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_spl3); -p_spl10 = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_spl10); -first(p_linear,5) +p_linear = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_linear); +p_spl3 = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_spl3); +p_spl10 = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_spl10); +first(p_linear, 5) # And plot them -pl = plot(evts.x,signal[1,1,:]) - lines!(p_linear.x,p_linear.yhat) - lines!(p_spl3.x,p_spl3.yhat) - lines!(p_spl10.x,p_spl10.yhat) +pl = plot(evts.x, signal[1, 1, :]) +lines!(p_linear.x, p_linear.yhat) +lines!(p_spl3.x, p_spl3.yhat) +lines!(p_spl10.x, p_spl10.yhat) pl # What we can clearly see here, that the linear effect (blue) underfits the data, the `spl(x,10)` overfits it, but the `spl(x,4)` fits it perfectly. @@ -72,46 +72,46 @@ typeof(term_spl) # a special type generated in Unfold.jl # # it has a field .fun which is the spline function. We can evaluate it at a point -const splFunction = Base.get_extension(Unfold,:UnfoldBSplineKitExt).splFunction -splFunction([0.2],term_spl) +const splFunction = Base.get_extension(Unfold, :UnfoldBSplineKitExt).splFunction +splFunction([0.2], term_spl) # each column of this 1 row matrix is a coefficient for our regression model. -lines(disallowmissing(splFunction([0.2],term_spl))[1,:]) +lines(disallowmissing(splFunction([0.2], term_spl))[1, :]) # Note: We have to use disallowmissing, because our splines return a `missing` whenever we ask it to return a value outside it's defined range, e.g.: -splFunction([-0.2],term_spl) +splFunction([-0.2], term_spl) # because it never has seen any data outside and can't extrapolate! # Okay back to our main issue: Let's plot the whole basis set -basisSet = splFunction(0.:0.01:1,term_spl) -basisSet = disallowmissing(basisSet[.!any(ismissing.(basisSet),dims=2)[:,1],:]) # remove missings -ax = Axis(Figure()[1,1]) -[lines!(ax,basisSet[:,k]) for k in 1:size(basisSet,2)] +basisSet = splFunction(0.0:0.01:1, term_spl) +basisSet = disallowmissing(basisSet[.!any(ismissing.(basisSet), dims = 2)[:, 1], :]) # remove missings +ax = Axis(Figure()[1, 1]) +[lines!(ax, basisSet[:, k]) for k = 1:size(basisSet, 2)] current_figure() # notice how we flipped the plot around, i.e. now on the x-axis we do not plot the coefficients, but the `x`-values # Now each line is one basis-function of the spline. # # Unfold returns us one coefficient per basis-function -β=coef(uf_spl10)[1,1,:] +β = coef(uf_spl10)[1, 1, :] β = Float64.(disallowmissing(β)) # but because we used an intercept, we have to do some remodelling in the basisSet -X = hcat(ones(size(basisSet,1)), basisSet[:,1:5], basisSet[:,7:end]) +X = hcat(ones(size(basisSet, 1)), basisSet[:, 1:5], basisSet[:, 7:end]) # now we can weight the spline by the basisfunction -weighted = (β .*X') +weighted = (β .* X') # Plotting them makes for a nice looking plot! -ax = Axis(Figure()[1,1]) -[lines!(weighted[k,:]) for k = 1:10] +ax = Axis(Figure()[1, 1]) +[lines!(weighted[k, :]) for k = 1:10] current_figure() # and sum them up -lines(sum(weighted,dims=1)[1,:]) -plot!(X*β,color="gray") #(same as matrixproduct X*β directly!) +lines(sum(weighted, dims = 1)[1, :]) +plot!(X * β, color = "gray") #(same as matrixproduct X*β directly!) current_figure() # And this is how you can think about splines diff --git a/docs/literate/explanations/window_length.jl b/docs/literate/explanations/window_length.jl index 06d811b5..9cfe4b84 100644 --- a/docs/literate/explanations/window_length.jl +++ b/docs/literate/explanations/window_length.jl @@ -23,65 +23,88 @@ set_theme!(theme_ggthemr(:fresh)) # # Init functions # First we need a function that simulates some continous data; conviently we can use UnfoldSim for this -function gen_data(rng,noiselevel,sfreq) - noise = PinkNoise(;noiselevel=noiselevel); - - dat,evts = UnfoldSim.predef_eeg(rng; +function gen_data(rng, noiselevel, sfreq) + noise = PinkNoise(; noiselevel = noiselevel) + + dat, evts = UnfoldSim.predef_eeg( + rng; sfreq = sfreq, - p1 = (p100(;sfreq=sfreq), @formula(0~1+condition),[5,0],Dict()), - n1 = (n170(;sfreq=sfreq), @formula(0~1+condition),[5,0],Dict()), - p3 = (p300(;sfreq=sfreq), @formula(0~1+continuous),[5,0],Dict()), - n_repeats=20,noise=noise) - return dat,evts + p1 = (p100(; sfreq = sfreq), @formula(0 ~ 1 + condition), [5, 0], Dict()), + n1 = (n170(; sfreq = sfreq), @formula(0 ~ 1 + condition), [5, 0], Dict()), + p3 = (p300(; sfreq = sfreq), @formula(0 ~ 1 + continuous), [5, 0], Dict()), + n_repeats = 20, + noise = noise, + ) + return dat, evts end; # Next a convience function to calculate the estimates -function calc_time_models(evts,dat,tWinList,sfreq) - mList = [] - for twindow = tWinList - m = fit(UnfoldModel,Dict(Any=>(@formula(0~1),firbasis(twindow,sfreq,string(twindow)))),evts,dat); - res = coeftable(m) - push!(mList,res) - end - return vcat(mList...) +function calc_time_models(evts, dat, tWinList, sfreq) + mList = [] + for twindow in tWinList + m = fit( + UnfoldModel, + Dict(Any => (@formula(0 ~ 1), firbasis(twindow, sfreq, string(twindow)))), + evts, + dat, + ) + res = coeftable(m) + push!(mList, res) + end + return vcat(mList...) end; # This is a convience function for plotting further down function basisToWin(basisname) - return parse.(Float64,[s[2] for s in split.(replace.(basisname,"("=>"",")"=>""),',')]) + return parse.( + Float64, + [s[2] for s in split.(replace.(basisname, "(" => "", ")" => ""), ',')], + ) end; # # Init variables -tWinList = [(-0.1,x) for x in [3,2.5,2,1.5,1,0.5]] +tWinList = [(-0.1, x) for x in [3, 2.5, 2, 1.5, 1, 0.5]] noiselevel = 8.5 sfreq = 250; # # Generate data and calculate estimates -dat,evts = gen_data(MersenneTwister(2),noiselevel,sfreq); +dat, evts = gen_data(MersenneTwister(2), noiselevel, sfreq); -res = calc_time_models(evts,dat,tWinList,sfreq); +res = calc_time_models(evts, dat, tWinList, sfreq); # We also append some additional information to the results dataframe res.tWin .= basisToWin(res.basisname); # For comparison lets also generate the ground truth of our data; this is a bit cumbersome and you don't have to care (too much) about it -dat_gt,evts_gt = UnfoldSim.predef_eeg(; - p1 = (p100(;sfreq=sfreq), @formula(0~1),[5],Dict()), - sfreq = sfreq, - n1 = (n170(;sfreq=sfreq), @formula(0~1),[5],Dict()), - p3 = (p300(;sfreq=sfreq), @formula(0~1),[5],Dict()), - n_repeats=1,noiselevel=0, return_epoched=true); -time_gt = range(0,length = length(dat_gt[:,1]),step=1/sfreq) +dat_gt, evts_gt = UnfoldSim.predef_eeg(; + p1 = (p100(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict()), + sfreq = sfreq, + n1 = (n170(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict()), + p3 = (p300(; sfreq = sfreq), @formula(0 ~ 1), [5], Dict()), + n_repeats = 1, + noiselevel = 0, + return_epoched = true, +); +time_gt = range(0, length = length(dat_gt[:, 1]), step = 1 / sfreq) df_gt = DataFrame( - basisname = reduce(vcat,fill.(unique(res.basisname), length(dat_gt[:,1]))), - channel = repeat([1], length(dat_gt[:,1])*length(unique(res.basisname))), - coefname = reduce(vcat,fill("GroundTruth", length(dat_gt[:,1])*length(unique(res.basisname)))), - estimate = repeat(dat_gt[:,1], length(unique(res.basisname))), - group = reduce(vcat,fill(nothing, length(dat_gt[:,1])*length(unique(res.basisname)))), - stderror = reduce(vcat,fill(nothing, length(dat_gt[:,1])*length(unique(res.basisname)))), + basisname = reduce(vcat, fill.(unique(res.basisname), length(dat_gt[:, 1]))), + channel = repeat([1], length(dat_gt[:, 1]) * length(unique(res.basisname))), + coefname = reduce( + vcat, + fill("GroundTruth", length(dat_gt[:, 1]) * length(unique(res.basisname))), + ), + estimate = repeat(dat_gt[:, 1], length(unique(res.basisname))), + group = reduce( + vcat, + fill(nothing, length(dat_gt[:, 1]) * length(unique(res.basisname))), + ), + stderror = reduce( + vcat, + fill(nothing, length(dat_gt[:, 1]) * length(unique(res.basisname))), + ), time = repeat(time_gt, length(unique(res.basisname))), - ); +); df_gt.tWin .= basisToWin(df_gt.basisname); @@ -93,32 +116,46 @@ res_gt = vcat(res, df_gt); # Plot figure f = Figure(); - + # Choose which data to plot -h_t = data(res) * mapping(:time,:estimate,color=:tWin,group=(:tWin,:coefname)=>(x,y)->string(x)*y); +h_t = + data(res) * mapping( + :time, + :estimate, + color = :tWin, + group = (:tWin, :coefname) => (x, y) -> string(x) * y, + ); # We use the following to plot some length indicator lines untWin = unique(res_gt.tWin) -segDF = DataFrame( - :x=>hcat(repeat([-0.1],length(untWin)),untWin)'[:], - :y=>repeat(reverse(1:length(untWin)),inner=2)) -segDF.tWin .= 0. +segDF = DataFrame( + :x => hcat(repeat([-0.1], length(untWin)), untWin)'[:], + :y => repeat(reverse(1:length(untWin)), inner = 2), +) +segDF.tWin .= 0.0 segDF.tWin[1:2:end] .= segDF.x[2:2:end] segDF.tWin[2:2:end] .= segDF.x[2:2:end] segDF.y = segDF.y .* 0.2 .+ 6; # Layer for indicator lines -h_l = data(@subset(segDF,:tWin .!=3.0))* mapping(:x,:y,color=:tWin,group=:tWin=>x->string.(x)); +h_l = + data(@subset(segDF, :tWin .!= 3.0)) * + mapping(:x, :y, color = :tWin, group = :tWin => x -> string.(x)); # Ground truth Layer -h_gt = data(df_gt) * mapping(:time,:estimate,group=(:tWin,:coefname)=>(x,y)->string(x)*y)*visual(Lines,linewidth=5,color=Colors.Gray(0.6)); +h_gt = + data(df_gt) * + mapping(:time, :estimate, group = (:tWin, :coefname) => (x, y) -> string(x) * y) * + visual(Lines, linewidth = 5, color = Colors.Gray(0.6)); # Add all visuals together and draw -h1 = h_gt + visual(Lines,colormap=get(ColorSchemes.Blues,0.3:0.01:1.2)) * (h_l + h_t) |> x->draw!(f[1,1],x,axis=(;xlabel="time [s]",ylabel="estimate [a.u.]")); +h1 = + h_gt + visual(Lines, colormap = get(ColorSchemes.Blues, 0.3:0.01:1.2)) * (h_l + h_t) |> + x -> draw!(f[1, 1], x, axis = (; xlabel = "time [s]", ylabel = "estimate [a.u.]")); # Add zero grid lines -h1 = hlines!(current_axis(),[0],color=Colors.Gray(0.8)); -h2 = vlines!(current_axis(),[0],color=Colors.Gray(0.8)); +h1 = hlines!(current_axis(), [0], color = Colors.Gray(0.8)); +h2 = vlines!(current_axis(), [0], color = Colors.Gray(0.8)); translate!(h1, 0, 0, -1); translate!(h2, 0, 0, -1); diff --git a/docs/make.jl b/docs/make.jl index 7aa31cee..5cb66a4a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,50 +7,56 @@ using Literate using Glob GENERATED = joinpath(@__DIR__, "src", "generated") -SOURCE = joinpath(@__DIR__,"literate") -for subfolder ∈ ["explanations","HowTo","tutorials"] - local SOURCE_FILES = Glob.glob(subfolder*"/*.jl", SOURCE) - foreach(fn -> Literate.markdown(fn, GENERATED*"/"*subfolder), SOURCE_FILES) +SOURCE = joinpath(@__DIR__, "literate") +for subfolder ∈ ["explanations", "HowTo", "tutorials"] + local SOURCE_FILES = Glob.glob(subfolder * "/*.jl", SOURCE) + foreach(fn -> Literate.markdown(fn, GENERATED * "/" * subfolder), SOURCE_FILES) end -makedocs(sitename="Unfold.jl Timeseries Analysis & Deconvolution", - #root = joinpath(dirname(pathof(Unfold)), "..", "docs"), - #prettyurls = get(ENV, "CI", nothing) == "true", - repo = Documenter.Remotes.GitHub("unfoldtoolbox", "Unfold.jl"), - pages = [ - "index.md", - "Installing Julia + Unfold.jl" => "installation.md", - "Tutorials"=>[ - - "Mass univariate LM" =>"tutorials/lm_mu.md", - "LM overlap correction" =>"tutorials/lm_overlap.md", - "Mass univariate Mixed Model" =>"tutorials/lmm_mu.md", - "LMM + overlap correction" =>"tutorials/lmm_overlap.md", - ], - "HowTo"=>[ - "Overlap: Multiple events"=>"HowTo/multiple_events.md", - "Import EEG with 🐍 PyMNE.jl"=>"HowTo/pymne.md" , - "Standard errors"=>"HowTo/standarderrors.md", - "Alternative Solvers (Robust, GPU, B2B)"=>"HowTo/custom_solvers.md", - "🐍 Calling Unfold.jl directly from Python" => "generated/HowTo/juliacall_unfold.md", - "P-values for mixedModels" => "HowTo/lmm_pvalues.md", - "Marginal effects (focus on non-linear predictors)" =>"generated/HowTo/effects.md", - #"Time domain basis functions"=>"generated/HowTo/timesplines.md", - "Save and load Unfold models" => "generated/HowTo/unfold_io.md", - ], - "Explanations"=>[ - - "About basisfunctions" => "./explanations/basisfunctions.md", - "Non-Linear effects" => "./generated/explanations/nonlinear_effects.md", - "Window Length Effect" => "./generated/explanations/window_length.md",], - "Reference"=>[ - "Overview of package extensions" => "references/extensions.md", - "Development environment" => "explanations/development.md", - "API: Types" => "references/types.md", - "API: Functions" => "references/functions.md"], - ]) - -deploydocs(; repo = "github.com/unfoldtoolbox/Unfold.jl", push_preview = true, devbranch = "main") +makedocs( + sitename = "Unfold.jl Timeseries Analysis & Deconvolution", + #root = joinpath(dirname(pathof(Unfold)), "..", "docs"), + #prettyurls = get(ENV, "CI", nothing) == "true", + repo = Documenter.Remotes.GitHub("unfoldtoolbox", "Unfold.jl"), + pages = [ + "index.md", + "Installing Julia + Unfold.jl" => "installation.md", + "Tutorials" => [ + "Mass univariate LM" => "tutorials/lm_mu.md", + "LM overlap correction" => "tutorials/lm_overlap.md", + "Mass univariate Mixed Model" => "tutorials/lmm_mu.md", + "LMM + overlap correction" => "tutorials/lmm_overlap.md", + ], + "HowTo" => [ + "Overlap: Multiple events" => "HowTo/multiple_events.md", + "Import EEG with 🐍 PyMNE.jl" => "HowTo/pymne.md", + "Standard errors" => "HowTo/standarderrors.md", + "Alternative Solvers (Robust, GPU, B2B)" => "HowTo/custom_solvers.md", + "🐍 Calling Unfold.jl directly from Python" => "generated/HowTo/juliacall_unfold.md", + "P-values for mixedModels" => "HowTo/lmm_pvalues.md", + "Marginal effects (focus on non-linear predictors)" => "generated/HowTo/effects.md", + #"Time domain basis functions"=>"generated/HowTo/timesplines.md", + "Save and load Unfold models" => "generated/HowTo/unfold_io.md", + ], + "Explanations" => [ + "About basisfunctions" => "./explanations/basisfunctions.md", + "Non-Linear effects" => "./generated/explanations/nonlinear_effects.md", + "Window Length Effect" => "./generated/explanations/window_length.md", + ], + "Reference" => [ + "Overview of package extensions" => "references/extensions.md", + "Development environment" => "explanations/development.md", + "API: Types" => "references/types.md", + "API: Functions" => "references/functions.md", + ], + ], +) + +deploydocs(; + repo = "github.com/unfoldtoolbox/Unfold.jl", + push_preview = true, + devbranch = "main", +) diff --git a/docs/nb_timesplines.jl b/docs/nb_timesplines.jl index 6d54e75e..28edb73c 100644 --- a/docs/nb_timesplines.jl +++ b/docs/nb_timesplines.jl @@ -6,47 +6,47 @@ using InteractiveUtils # ╔═╡ 25d21156-d006-11eb-3655-61be76e7db93 begin -import Pkg; - Pkg.activate("../../../") + import Pkg + Pkg.activate("../../../") end # ╔═╡ 34aefa9a-e034-4293-b53d-cc1dc349ecc7 let - using Revise - using Unfold - using DataFramesMeta + using Revise + using Unfold + using DataFramesMeta end # ╔═╡ c73bf3ed-4af3-4e5a-997e-6399c9b85bbe -using StatsModels,PlutoUI,DataFrames,StatsPlots +using StatsModels, PlutoUI, DataFrames, StatsPlots # ╔═╡ d1e1de70-7e37-46a4-8ce2-3ca591bb84cf -include(pathof(Unfold)*"../../../test/test_utilities.jl") +include(pathof(Unfold) * "../../../test/test_utilities.jl") # ╔═╡ 2e9ec6e0-03c6-41f6-bd3b-a44eaed5b553 -data,evts = loadtestdata("test_case_4a"); +data, evts = loadtestdata("test_case_4a"); # ╔═╡ 0b27e0bb-7aee-4672-8501-6096539c7ad7 -evts.continuousA = rand(size(evts,1)) +evts.continuousA = rand(size(evts, 1)) # ╔═╡ 97443eb8-ac67-4f8e-9e8a-eb09176c8516 # ╔═╡ 2a5d5493-b88d-473f-8d81-9e9240b3ffe0 #f = @formula 0~1+conditionA+continuousA # 1 -f = @formula 0~1+continuousA # 1 +f = @formula 0 ~ 1 + continuousA # 1 # ╔═╡ 5d23b74c-988d-48bc-9c30-c2af0dbabdb4 # ╔═╡ 3f2b7a92-a064-4be6-8a39-cabaa2453777 -basisfunction = Unfold.splinebasis(τ=(-1,1),sfreq=20,nsplines=10,name="basisA") +basisfunction = Unfold.splinebasis(τ = (-1, 1), sfreq = 20, nsplines = 10, name = "basisA") # ╔═╡ 8f32dbf3-8958-4e67-8e25-5f24d48058e1 Unfold.splinebasis # ╔═╡ 7ef24e5a-63df-4e99-a021-7119d7e8d3bf -fir = firbasis(τ=(-1,1),sfreq=20,name="basisB") +fir = firbasis(τ = (-1, 1), sfreq = 20, name = "basisB") # ╔═╡ 12c617c1-1994-4efe-b8bc-23fd2325c2ca fir.kernel(2) @@ -67,19 +67,25 @@ size(basisfunction.colnames) evts # ╔═╡ c4227637-81b4-4ea7-a7b7-741860fef8c3 -m = fit(UnfoldModel,Dict("eventA"=>(f,fir),"eventB"=>(f,basisfunction)),evts,data,eventcolumn="type"); +m = fit( + UnfoldModel, + Dict("eventA" => (f, fir), "eventB" => (f, basisfunction)), + evts, + data, + eventcolumn = "type", +); # ╔═╡ cd19784c-686e-41c5-bd32-cdf793cecf03 res = coeftable(m) # ╔═╡ 29e875d3-45b9-4a6a-aaa6-c6eaf35244c5 -@df res plot(:time,:estimate,group=(:basisname,:coefname)) +@df res plot(:time, :estimate, group = (:basisname, :coefname)) # ╔═╡ b669f55b-9121-4d66-9c82-f89bbfe0b2df Unfold.coef(m) # ╔═╡ 64cf7e78-2677-4551-9491-82c9f560ed61 - res +res # ╔═╡ 02b8ff44-c514-492e-bf94-0cbf44431097 x = [1] @@ -89,70 +95,74 @@ x = [1] # ╔═╡ 1ee404e1-3039-410d-99de-ce713e751280 let - basisnames = Unfold.get_basis_name(m) - betaOut = [] - for (k,b) = enumerate(unique(basisnames)) - - ix = basisnames .== b - - betaOut = vcat(betaOut,designmatrix(m).formulas[k].rhs.basisfunction.kernel(0) * m.modelfit.estimate[:,ix]') - end - plot(Float64.(betaOut)) + basisnames = Unfold.get_basis_name(m) + betaOut = [] + for (k, b) in enumerate(unique(basisnames)) + + ix = basisnames .== b + + betaOut = vcat( + betaOut, + designmatrix(m).formulas[k].rhs.basisfunction.kernel(0) * + m.modelfit.estimate[:, ix]', + ) + end + plot(Float64.(betaOut)) end # ╔═╡ 5466ec0b-d28a-4be7-b686-cb5ad7fea92b let - basisnames = Unfold.get_basis_name(m) - k = 2 - b = unique(basisnames)[k] - - - bf = designmatrix(m).formulas[k].rhs.basisfunction - ix = basisnames .== b - - hcat(bf.kernel.([0,0]) ...) * m.modelfit.estimate[:,ix]' + basisnames = Unfold.get_basis_name(m) + k = 2 + b = unique(basisnames)[k] + + + bf = designmatrix(m).formulas[k].rhs.basisfunction + ix = basisnames .== b + + hcat(bf.kernel.([0, 0])...) * m.modelfit.estimate[:, ix]' end # ╔═╡ 033cf2e0-a3ca-40e3-8f27-e88e017ca8de size(Unfold.times(basisfunction)) # ╔═╡ a8238bd7-0198-45d2-89bd-d1434c537b2a - designmatrix(m).formulas[2].rhs.basisfunction.kernel(0) +designmatrix(m).formulas[2].rhs.basisfunction.kernel(0) # ╔═╡ fe5a826f-aede-4e91-b6cd-ee7fca33ce31 zeros() # ╔═╡ 8b9d8e32-d2ca-49db-986f-851ffabc3a3e -function undoBasis(d,m) - bname = unique(d.basisname) # the current coefficient basename - @assert(length(bname)==1) - @assert(length(unique(d.channel))==1) - bnames = unique(Unfold.get_basis_name(m)) # all model basenames - - # find out which formula / basisfunction we have - k = findfirst(bname .== bnames) - bf = designmatrix(m).formulas[k].rhs.basisfunction - - - estimate= bf.kernel.(0) *d.estimate - @show size(bf.kernel.(0)) - @show size(Unfold.times(bf)) - return DataFrame(:time=>Unfold.times(bf),:estimate=>estimate) +function undoBasis(d, m) + bname = unique(d.basisname) # the current coefficient basename + @assert(length(bname) == 1) + @assert(length(unique(d.channel)) == 1) + bnames = unique(Unfold.get_basis_name(m)) # all model basenames + + # find out which formula / basisfunction we have + k = findfirst(bname .== bnames) + bf = designmatrix(m).formulas[k].rhs.basisfunction + + + estimate = bf.kernel.(0) * d.estimate + @show size(bf.kernel.(0)) + @show size(Unfold.times(bf)) + return DataFrame(:time => Unfold.times(bf), :estimate => estimate) end # ╔═╡ 350a7f4b-f019-4bf6-8878-2cbb0cf82005 begin - gd = groupby(res,[:basisname,:channel,:group,:coefname]) - a = combine(gd) do d - return undoBasis(d,m) - end - a + gd = groupby(res, [:basisname, :channel, :group, :coefname]) + a = combine(gd) do d + return undoBasis(d, m) + end + a end # ╔═╡ a876f554-a034-461e-9522-490ca4bab5f8 -@df a plot(:time,:estimate,group=(:basisname,:coefname)) +@df a plot(:time, :estimate, group = (:basisname, :coefname)) # ╔═╡ 1d5b8b05-cdc6-4058-9bc1-bd52e6e3b444 Unfold.get_basis_name(m) diff --git a/docs/nb_timesplines_dev.jl b/docs/nb_timesplines_dev.jl index 39379804..47e92df5 100644 --- a/docs/nb_timesplines_dev.jl +++ b/docs/nb_timesplines_dev.jl @@ -7,7 +7,14 @@ using InteractiveUtils # This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). macro bind(def, element) quote - local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local iv = try + Base.loaded_modules[Base.PkgId( + Base.UUID("6e696c72-6542-2067-7265-42206c756150"), + "AbstractPlutoDingetjes", + )].Bonds.initial_value + catch + b -> missing + end local el = $(esc(element)) global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) el @@ -21,7 +28,7 @@ import Pkg; Pkg.activate("/store/users/ehinger/unfoldjl_dev/") # ╔═╡ cf3dcd54-cfa0-11eb-2e99-756f9df71a4d -using Revise,Unfold,Plots,BSplines +using Revise, Unfold, Plots, BSplines # ╔═╡ a9e6971f-8a50-4564-ace5-1fabe65b1522 import PlutoUI @@ -36,50 +43,50 @@ import PlutoUI # ╔═╡ 89eace50-c0f5-4951-821f-cf010cdb9726 -@bind nsplines PlutoUI.Slider(2:10; default=8, show_value=true) +@bind nsplines PlutoUI.Slider(2:10; default = 8, show_value = true) # ╔═╡ cb62fb85-364d-4142-b28b-c8b0871b4ccc -@bind srate PlutoUI.Slider(1:50; default=10, show_value=true) +@bind srate PlutoUI.Slider(1:50; default = 10, show_value = true) # ╔═╡ 8496a69b-20bc-45c2-92cb-179d0ad7127c -@bind τ1 PlutoUI.Slider(-5:0.1:2; default=-.5, show_value=true) +@bind τ1 PlutoUI.Slider(-5:0.1:2; default = -0.5, show_value = true) # ╔═╡ c8285ed7-bffb-4047-89ff-22a664c19a78 -@bind τ2 PlutoUI.Slider(-2:0.1:5; default=2, show_value=true) +@bind τ2 PlutoUI.Slider(-2:0.1:5; default = 2, show_value = true) # ╔═╡ 9df0d50d-2a1a-4431-9725-12178c824d62 -τ = [τ1,τ2] +τ = [τ1, τ2] # ╔═╡ 94a258ea-febd-4793-89c4-d3755e82c48c begin - function splinebasis(τ,sfreq,nsplines,name::String) - τ = Unfold.round_times(τ,sfreq) - times =range(τ[1],stop=τ[2],step=1 ./sfreq) - kernel=e->splinekernel(e,times,nsplines) - type = "splinebasis" - - shiftOnset = Int64(floor(τ[1] * sfreq)) - - return Unfold.BasisFunction(kernel,times,times,type,name,shiftOnset) - end - - - - function spl_breakpoints(times,nsplines) - # calculate the breakpoints, evenly spaced - return collect(range(minimum(times),stop=maximum(times),length=nsplines)) - end - function splinekernel(e,times,nsplines) - breakpoints = spl_breakpoints(times,nsplines) - basis = BSplineBasis(4,breakpoints) # 4= cubic - return Unfold.splFunction(times,basis) - end - spl=splinebasis(τ,srate,nsplines,"test") - + function splinebasis(τ, sfreq, nsplines, name::String) + τ = Unfold.round_times(τ, sfreq) + times = range(τ[1], stop = τ[2], step = 1 ./ sfreq) + kernel = e -> splinekernel(e, times, nsplines) + type = "splinebasis" + + shiftOnset = Int64(floor(τ[1] * sfreq)) + + return Unfold.BasisFunction(kernel, times, times, type, name, shiftOnset) + end + + + + function spl_breakpoints(times, nsplines) + # calculate the breakpoints, evenly spaced + return collect(range(minimum(times), stop = maximum(times), length = nsplines)) + end + function splinekernel(e, times, nsplines) + breakpoints = spl_breakpoints(times, nsplines) + basis = BSplineBasis(4, breakpoints) # 4= cubic + return Unfold.splFunction(times, basis) + end + spl = splinebasis(τ, srate, nsplines, "test") + end # ╔═╡ 3f612e26-aba8-46ed-af19-67cd1be672fa -plot(spl.times,spl.kernel(0)) +plot(spl.times, spl.kernel(0)) # ╔═╡ d375b64c-0e64-4b39-8c3f-5a14d365877d diff --git a/docs/nb_unfold-overlap.jl b/docs/nb_unfold-overlap.jl index a0320056..aeffe368 100644 --- a/docs/nb_unfold-overlap.jl +++ b/docs/nb_unfold-overlap.jl @@ -15,16 +15,16 @@ end # ╔═╡ 65edde2e-d03e-11eb-300b-897bdc66d875 begin - using Unfold - using PlutoUI - using StatsPlots - using DataFrames - using StatsModels + using Unfold + using PlutoUI + using StatsPlots + using DataFrames + using StatsModels end # ╔═╡ 0a821c04-1aff-4fca-aab8-07a4d450838c begin - using Random + using Random end # ╔═╡ bec07615-8ca3-4251-9e44-bf5044828274 @@ -42,10 +42,10 @@ Setup might take a while because the libraries are very large. It might be bette md"""### Loading Data & Setting up Unfold""" # ╔═╡ 76afd758-ac50-4f3c-a484-7f268ab059ef -data,evts = loadtestdata("test_case_3b"); +data, evts = loadtestdata("test_case_3b"); # ╔═╡ c7fb216e-d177-4623-a135-db4cd856d4ad -f = @formula 0~1+conditionA+continuousA +f = @formula 0 ~ 1 + conditionA + continuousA # ╔═╡ a55058c9-3b56-447b-8f1a-f7319a29f478 md"""### Fitting the models""" @@ -54,24 +54,24 @@ md"""### Fitting the models""" md"""### Plotting the results""" # ╔═╡ 9a4c715a-7a51-411d-a354-f79ec2369cb7 -let - md"""change window size τ = (-0.3, $(@bind τ2 Slider(0:0.1:5,default=2, show_value=true))) - """ +let + md"""change window size τ = (-0.3, $(@bind τ2 Slider(0:0.1:5,default=2, show_value=true))) + """ end # ╔═╡ ba2f0c48-f5f8-40eb-a45b-441b8134ffb1 -τ = (-0.3,τ2) +τ = (-0.3, τ2) # ╔═╡ f99efde3-bd34-46e2-ad1c-1a97ad866b79 -fir = firbasis(τ=τ,sfreq=20,name="basisA"); +fir = firbasis(τ = τ, sfreq = 20, name = "basisA"); # ╔═╡ dcbdef13-675e-47f7-b484-d340f0e55934 -f_dict = Dict("eventA"=>(f,fir)) # combine eventname, formula & basis +f_dict = Dict("eventA" => (f, fir)) # combine eventname, formula & basis # ╔═╡ ed61f099-343c-47d4-9af4-c4486c93fd4a -let - md"""change noise: σ = $(@bind σ Slider(0:0.1:5,default=0, show_value=true)) - """ +let + md"""change noise: σ = $(@bind σ Slider(0:0.1:5,default=0, show_value=true)) + """ end # ╔═╡ 09705728-c6e5-41a5-9859-2c7e28828d5b @@ -80,28 +80,44 @@ data_noise = data .+ σ .* randn(size(data)); # ╔═╡ b31a80a5-cca2-4c69-82f5-2b477f5b62be begin - plot(data_noise[1:600],title="Simulated Data with overlap") - vline!(evts.latency[1:20],legend=false) + plot(data_noise[1:600], title = "Simulated Data with overlap") + vline!(evts.latency[1:20], legend = false) end # ╔═╡ 3f4dd059-e275-483d-af42-a09a2795fca4 -data_e,times = Unfold.epoch(data=data_noise,tbl=evts,τ=τ,sfreq=10); +data_e, times = Unfold.epoch(data = data_noise, tbl = evts, τ = τ, sfreq = 10); # ╔═╡ 9aa42667-2629-48e1-b7fe-979dabc28f96 # non-overlapping (mass univariate) -m_e,res_e = fit(UnfoldLinearModel,f,evts,data_e,times); +m_e, res_e = fit(UnfoldLinearModel, f, evts, data_e, times); # ╔═╡ bce42119-33d7-4e10-ac5e-b37028381b09 # overlap-corrected -m,res = fit(UnfoldLinearModel,f_dict,evts,data_noise,eventcolumn="type"); +m, res = fit(UnfoldLinearModel, f_dict, evts, data_noise, eventcolumn = "type"); # ╔═╡ de267143-c029-4103-9fd7-ada153588ab6 begin - - -p1 = @df res_e plot(:colname_basis,:estimate,group=:coefname,title="Without ...",ylims=(-3,5),legend=false,xlims=(-0.3,3)) - p2 = @df res plot(:colname_basis,:estimate,group=:coefname,title="... with overlap correction",ylims=(-3,5),xlims=(-0.3,3),legend=:bottomright) - plot(p1,p2,layout=(1,2)) + + + p1 = @df res_e plot( + :colname_basis, + :estimate, + group = :coefname, + title = "Without ...", + ylims = (-3, 5), + legend = false, + xlims = (-0.3, 3), + ) + p2 = @df res plot( + :colname_basis, + :estimate, + group = :coefname, + title = "... with overlap correction", + ylims = (-3, 5), + xlims = (-0.3, 3), + legend = :bottomright, + ) + plot(p1, p2, layout = (1, 2)) end # ╔═╡ Cell order: diff --git a/docs/nb_unfold-overlap_extended.jl b/docs/nb_unfold-overlap_extended.jl index 0710c76a..567e6114 100644 --- a/docs/nb_unfold-overlap_extended.jl +++ b/docs/nb_unfold-overlap_extended.jl @@ -15,22 +15,22 @@ end # ╔═╡ 65edde2e-d03e-11eb-300b-897bdc66d875 begin -import Pkg; -Pkg.activate(mktempdir()) - Pkg.add(url="https://github.com/unfoldtoolbox/Unfold.jl") - Pkg.add("PlutoUI") - Pkg.add("StatsPlots") - Pkg.add("DataFrames") - Pkg.add("StatsModels") - Pkg.add("DSP") + import Pkg + Pkg.activate(mktempdir()) + Pkg.add(url = "https://github.com/unfoldtoolbox/Unfold.jl") + Pkg.add("PlutoUI") + Pkg.add("StatsPlots") + Pkg.add("DataFrames") + Pkg.add("StatsModels") + Pkg.add("DSP") end # ╔═╡ 0a821c04-1aff-4fca-aab8-07a4d450838c let - using Revise - using Unfold - using StatsModels,PlutoUI,StatsPlots,DataFrames,Random - using DSP + using Revise + using Unfold + using StatsModels, PlutoUI, StatsPlots, DataFrames, Random + using DSP end # ╔═╡ bec07615-8ca3-4251-9e44-bf5044828274 @@ -48,10 +48,10 @@ Setup might take a while because the libraries are very large. It might be bette md"""### Loading Data & Setting up Unfold""" # ╔═╡ 76afd758-ac50-4f3c-a484-7f268ab059ef -data,evts = loadtestdata("test_case_3b"); +data, evts = loadtestdata("test_case_3b"); # ╔═╡ c7fb216e-d177-4623-a135-db4cd856d4ad -f = @formula 0~1+conditionA+continuousA +f = @formula 0 ~ 1 + conditionA + continuousA # ╔═╡ a55058c9-3b56-447b-8f1a-f7319a29f478 md"""### Fitting the models""" @@ -60,66 +60,82 @@ md"""### Fitting the models""" md"""### Plotting the results""" # ╔═╡ 9a4c715a-7a51-411d-a354-f79ec2369cb7 -let - md"""change window size τ = (-0.3, $(@bind τ2 Slider(0:0.1:5,default=2, show_value=true))) - """ +let + md"""change window size τ = (-0.3, $(@bind τ2 Slider(0:0.1:5,default=2, show_value=true))) + """ end # ╔═╡ ba2f0c48-f5f8-40eb-a45b-441b8134ffb1 -τ = (-0.3,τ2) +τ = (-0.3, τ2) # ╔═╡ f99efde3-bd34-46e2-ad1c-1a97ad866b79 -fir_basis = firbasis(τ=τ,sfreq=20,name="basisA"); +fir_basis = firbasis(τ = τ, sfreq = 20, name = "basisA"); # ╔═╡ dcbdef13-675e-47f7-b484-d340f0e55934 -f_dict = Dict("eventA"=>(f,fir_basis)) # combine eventname, formula & basis +f_dict = Dict("eventA" => (f, fir_basis)) # combine eventname, formula & basis # ╔═╡ ed61f099-343c-47d4-9af4-c4486c93fd4a -let - md"""change noise: σ = $(@bind σ Slider(0:0.1:5,default=0, show_value=true)) - """ +let + md"""change noise: σ = $(@bind σ Slider(0:0.1:5,default=0, show_value=true)) + """ end # ╔═╡ 09705728-c6e5-41a5-9859-2c7e28828d5b # add some noise #data_noise = data .+ σ .* randn(size(data)); -data_noise = data .+ σ .* rand(PinkGaussian(size(data,1))); +data_noise = data .+ σ .* rand(PinkGaussian(size(data, 1))); # ╔═╡ 1f3419f0-02e3-4346-8ac4-026cb76fa0fe - md"""filter: $(@bind low_cutoff Slider(0.01:0.01:2,default=0.01, show_value=true)) - """ +md"""filter: $(@bind low_cutoff Slider(0.01:0.01:2,default=0.01, show_value=true)) +""" # ╔═╡ 3e6f413d-f0b3-434f-a2bd-e002592247f8 -hpf = fir(501,low_cutoff; fs=20); +hpf = fir(501, low_cutoff; fs = 20); # ╔═╡ b2893e0a-5712-4af2-923d-e80078f7277c data_filter = filtfilt(hpf, data_noise); # ╔═╡ b31a80a5-cca2-4c69-82f5-2b477f5b62be begin - plot(data_filter[1:400],title="Simulated Data with overlap") - plot!(data_noise[1:400],title="Simulated Data with overlap") - vline!(evts.latency[1:10],legend=false) + plot(data_filter[1:400], title = "Simulated Data with overlap") + plot!(data_noise[1:400], title = "Simulated Data with overlap") + vline!(evts.latency[1:10], legend = false) end # ╔═╡ 3f4dd059-e275-483d-af42-a09a2795fca4 -data_e,times = Unfold.epoch(data=data_filter,tbl=evts,τ=τ,sfreq=20); +data_e, times = Unfold.epoch(data = data_filter, tbl = evts, τ = τ, sfreq = 20); # ╔═╡ 9aa42667-2629-48e1-b7fe-979dabc28f96 # non-overlapping (mass univariate) -m_e,res_e = fit(UnfoldLinearModel,f,evts,data_e,times); +m_e, res_e = fit(UnfoldLinearModel, f, evts, data_e, times); # ╔═╡ bce42119-33d7-4e10-ac5e-b37028381b09 # overlap-corrected -m,res = fit(UnfoldLinearModel,f_dict,evts,data_filter,eventcolumn="type"); +m, res = fit(UnfoldLinearModel, f_dict, evts, data_filter, eventcolumn = "type"); # ╔═╡ de267143-c029-4103-9fd7-ada153588ab6 begin - - -p1 = @df res_e plot(:colname_basis,:estimate,group=:coefname,title="Without ...",ylims=(-3,5),legend=false,xlims=(-0.3,5)) - p2 = @df res plot(:colname_basis,:estimate,group=:coefname,title="... with overlap correction",ylims=(-3,5),xlims=(-0.3,5),legend=:bottomright) - plot(p1,p2,layout=(1,2)) + + + p1 = @df res_e plot( + :colname_basis, + :estimate, + group = :coefname, + title = "Without ...", + ylims = (-3, 5), + legend = false, + xlims = (-0.3, 5), + ) + p2 = @df res plot( + :colname_basis, + :estimate, + group = :coefname, + title = "... with overlap correction", + ylims = (-3, 5), + xlims = (-0.3, 5), + legend = :bottomright, + ) + plot(p1, p2, layout = (1, 2)) end # ╔═╡ Cell order: diff --git a/docs/notebooks/nb_nonlinear.jl b/docs/notebooks/nb_nonlinear.jl index 0b7a83d1..0f4fb46f 100644 --- a/docs/notebooks/nb_nonlinear.jl +++ b/docs/notebooks/nb_nonlinear.jl @@ -6,57 +6,57 @@ using InteractiveUtils # ╔═╡ 1eec25bc-8040-11ec-024c-c54939c335ac begin - - using CairoMakie - using DataFrames - using Random -using Unfold -using Colors + + using CairoMakie + using DataFrames + using Random + using Unfold + using Colors end # ╔═╡ f926f2a9-d188-4c65-80c7-0e500d9c3b9a begin - rng = MersenneTwister(2) -n = 20 -evts = DataFrame(:x=>rand(rng,n)) -signal = -(3*(evts.x .-0.5)).^2 .+ 0.5 .* rand(rng,n) - signal = reshape(signal,length(signal),1,1) - signal = permutedims(signal,[3,2,1]) + rng = MersenneTwister(2) + n = 20 + evts = DataFrame(:x => rand(rng, n)) + signal = -(3 * (evts.x .- 0.5)) .^ 2 .+ 0.5 .* rand(rng, n) + signal = reshape(signal, length(signal), 1, 1) + signal = permutedims(signal, [3, 2, 1]) end; # ╔═╡ 248e07ea-bf5b-4c88-97c6-e29d0576bdaa begin - design_spl3 = Dict(Any=> (@formula(0~1+spl(x,3)),[0])) - uf_spl3 = fit(UnfoldModel,design_spl3,evts,signal); + design_spl3 = Dict(Any => (@formula(0 ~ 1 + spl(x, 3)), [0])) + uf_spl3 = fit(UnfoldModel, design_spl3, evts, signal) end; # ╔═╡ fea2e408-baa5-426c-a125-6ef68953abce begin - design_lin = Dict(Any=> (@formula(0~1+x),[0])) - uf_lin = fit(UnfoldModel,design_lin,evts,signal); + design_lin = Dict(Any => (@formula(0 ~ 1 + x), [0])) + uf_lin = fit(UnfoldModel, design_lin, evts, signal) end; # ╔═╡ f30a9476-b884-44e4-b2e9-b6ca722e52da begin - design_spl10 = Dict(Any=> (@formula(0~1+spl(x,10)),[0])) - uf_spl10 = fit(UnfoldModel,design_spl10,evts,signal); + design_spl10 = Dict(Any => (@formula(0 ~ 1 + spl(x, 10)), [0])) + uf_spl10 = fit(UnfoldModel, design_spl10, evts, signal) end; # ╔═╡ 144b86fc-b860-4f12-9250-0d7926e68cc5 begin - p_3 = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_spl3); - p_10 = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_spl10); - p_l = Unfold.effects(Dict(:x => range(0,stop=1,length=100)),uf_lin); + p_3 = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_spl3) + p_10 = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_spl10) + p_l = Unfold.effects(Dict(:x => range(0, stop = 1, length = 100)), uf_lin) end; # ╔═╡ 8633dc2f-cde4-45ca-8378-d4703804d02f begin - pl = plot(evts.x,signal[1,1,:]) - lines!(p_l.x,p_l.yhat) - lines!(p_3.x,p_3.yhat) - lines!(p_10.x,p_10.yhat) + pl = plot(evts.x, signal[1, 1, :]) + lines!(p_l.x, p_l.yhat) + lines!(p_3.x, p_3.yhat) + lines!(p_10.x, p_10.yhat) - pl + pl end # ╔═╡ 71e73918-f6c5-4f48-adf2-74875c32833c diff --git a/docs/run_liveserver.jl b/docs/run_liveserver.jl index efb6d45b..a7f5433f 100644 --- a/docs/run_liveserver.jl +++ b/docs/run_liveserver.jl @@ -1,7 +1,10 @@ using LiveServer -dr = filter(isdir,readdir(joinpath("src","generated"),join=true)) -push!(dr,"./build") -servedocs(skip_dirs=dr,literate_dir=joinpath("literate"),foldername=".",host="0.0.0.0") - - +dr = filter(isdir, readdir(joinpath("src", "generated"), join = true)) +push!(dr, "./build") +servedocs( + skip_dirs = dr, + literate_dir = joinpath("literate"), + foldername = ".", + host = "0.0.0.0", +) diff --git a/ext/UnfoldBSplineKitExt/UnfoldBSplineKitExt.jl b/ext/UnfoldBSplineKitExt/UnfoldBSplineKitExt.jl index 59229c24..abb673d1 100644 --- a/ext/UnfoldBSplineKitExt/UnfoldBSplineKitExt.jl +++ b/ext/UnfoldBSplineKitExt/UnfoldBSplineKitExt.jl @@ -1,8 +1,8 @@ module UnfoldBSplineKitExt using Unfold - using BSplineKit - using StatsModels -import StatsModels: termvars,width,coefnames,modelcols,apply_schema +using BSplineKit +using StatsModels +import StatsModels: termvars, width, coefnames, modelcols, apply_schema import Base: show using Statistics using DataFrames @@ -10,15 +10,16 @@ using Effects using SparseArrays import Unfold: AbstractSplineTerm - include("basisfunctions.jl") - include("splinepredictors.jl") +include("basisfunctions.jl") +include("splinepredictors.jl") ## Effects -Effects._trmequal(t1::AbstractSplineTerm,t2::AbstractTerm) = Effects._symequal(t1.term,t2) -Effects._trmequal(t1::AbstractSplineTerm,t2::AbstractSplineTerm) = Effects._symequal(t1.term,t2.term) +Effects._trmequal(t1::AbstractSplineTerm, t2::AbstractTerm) = Effects._symequal(t1.term, t2) +Effects._trmequal(t1::AbstractSplineTerm, t2::AbstractSplineTerm) = + Effects._symequal(t1.term, t2.term) -Effects._trmequal(t1::AbstractTerm,t2::AbstractSplineTerm) = Effects._symequal(t1,t2.term) -Effects._symequal(t1::AbstractTerm,t2::AbstractSplineTerm) = Effects._symequal(t1,t2.term) +Effects._trmequal(t1::AbstractTerm, t2::AbstractSplineTerm) = Effects._symequal(t1, t2.term) +Effects._symequal(t1::AbstractTerm, t2::AbstractSplineTerm) = Effects._symequal(t1, t2.term) -end \ No newline at end of file +end diff --git a/ext/UnfoldBSplineKitExt/basisfunctions.jl b/ext/UnfoldBSplineKitExt/basisfunctions.jl index bfa38398..14520dcc 100644 --- a/ext/UnfoldBSplineKitExt/basisfunctions.jl +++ b/ext/UnfoldBSplineKitExt/basisfunctions.jl @@ -1,6 +1,7 @@ splinebasis(; τ, sfreq, nsplines, name) = Unfold.splinebasis(τ, sfreq, nsplines, name) -splinebasis(τ, sfreq, nsplines) = Unfold.splinebasis(τ, sfreq, nsplines, "basis_" * string(rand(1:10000))) +splinebasis(τ, sfreq, nsplines) = + Unfold.splinebasis(τ, sfreq, nsplines, "basis_" * string(rand(1:10000))) @@ -26,4 +27,4 @@ function splinekernel(e, times, nsplines) breakpoints = spl_breakpoints(times, nsplines) basis = BSplineKit.BSplineBasis(BSplineOrder(4), breakpoints) # 4= cubic return sparse(splFunction(times, basis)) -end \ No newline at end of file +end diff --git a/ext/UnfoldBSplineKitExt/splinepredictors.jl b/ext/UnfoldBSplineKitExt/splinepredictors.jl index c622eb13..1c7ee6a7 100644 --- a/ext/UnfoldBSplineKitExt/splinepredictors.jl +++ b/ext/UnfoldBSplineKitExt/splinepredictors.jl @@ -24,7 +24,7 @@ Note: that due to the boundary condition (`natural`) spline, we repeat the bound """ function genSpl_breakpoints(p::AbstractSplineTerm, x) - p = range(0.0, length=p.df - 2, stop=1.0) + p = range(0.0, length = p.df - 2, stop = 1.0) breakpoints = quantile(x, p) return breakpoints end @@ -34,7 +34,7 @@ In the circular case, we do not use quantiles, (circular quantiles are difficult """ function genSpl_breakpoints(p::PeriodicBSplineTerm, x) # periodic case - - return range(p.low, p.high, length=p.df + 2) + return range(p.low, p.high, length = p.df + 2) end """ @@ -63,7 +63,9 @@ function spl_fillMat!(bs::BSplineBasis, large::Matrix, x::AbstractVector) ix = x .< bnds[1] .|| x .> bnds[2] if sum(ix) != 0 - @warn("spline prediction outside of possible range putting those values to missing.\n `findfirst(Out-Of-Bound-value)` is x=$(x[findfirst(ix)]), with bounds: $bnds") + @warn( + "spline prediction outside of possible range putting those values to missing.\n `findfirst(Out-Of-Bound-value)` is x=$(x[findfirst(ix)]), with bounds: $bnds" + ) large[ix, :] .= missing end @@ -98,17 +100,18 @@ Unfold.spl(x, df) = 0 # fallback # make a nice call if the function is called via REPL Unfold.spl(t::Symbol, d::Int) = BSplineTerm(term(t), d, 4, []) -Unfold.circspl(t::Symbol, d::Int, low, high) = PeriodicBSplineTerm(term(t), term(d), 4, low, high) +Unfold.circspl(t::Symbol, d::Int, low, high) = + PeriodicBSplineTerm(term(t), term(d), 4, low, high) """ Construct a BSplineTerm, if breakpoints/basis are not defined yet, put to `nothing` """ -function BSplineTerm(term, df, order=4) +function BSplineTerm(term, df, order = 4) @assert df > 3 "Minimal degrees of freedom has to be 4" BSplineTerm(term, df, order, []) end -function BSplineTerm(term, df::ConstantTerm, order=4) +function BSplineTerm(term, df::ConstantTerm, order = 4) BSplineTerm(term, df.n, order, []) end @@ -116,7 +119,14 @@ end function PeriodicBSplineTerm(term, df, low, high) PeriodicBSplineTerm(term, df, 4, low, high) end -function PeriodicBSplineTerm(term::AbstractTerm, df::ConstantTerm, order, low::ConstantTerm, high::ConstantTerm, breakvec) +function PeriodicBSplineTerm( + term::AbstractTerm, + df::ConstantTerm, + order, + low::ConstantTerm, + high::ConstantTerm, + breakvec, +) PeriodicBSplineTerm(term, df.n, order, low.n, high.n, breakvec) end function PeriodicBSplineTerm(term, df, order, low, high) @@ -124,7 +134,8 @@ function PeriodicBSplineTerm(term, df, order, low, high) end Base.show(io::IO, p::BSplineTerm) = print(io, "spl($(p.term), $(p.df))") -Base.show(io::IO, p::PeriodicBSplineTerm) = print(io, "circspl($(p.term), $(p.df),$(p.low):$(p.high))") +Base.show(io::IO, p::PeriodicBSplineTerm) = + print(io, "circspl($(p.term), $(p.df),$(p.low):$(p.high))") function StatsModels.apply_schema( t::FunctionTerm{typeof(Unfold.spl)}, @@ -175,7 +186,8 @@ function StatsModels.apply_schema( return construct_spline(t, term) end construct_spline(t::BSplineTerm, term) = BSplineTerm(term, t.df, t.order) -construct_spline(t::PeriodicBSplineTerm, term) = PeriodicBSplineTerm(term, t.df, t.order, t.low, t.high) +construct_spline(t::PeriodicBSplineTerm, term) = + PeriodicBSplineTerm(term, t.df, t.order, t.low, t.high) function StatsModels.modelcols(p::AbstractSplineTerm, d::NamedTuple) diff --git a/ext/UnfoldKrylovExt.jl b/ext/UnfoldKrylovExt.jl index 7acb5d97..1d6f48aa 100644 --- a/ext/UnfoldKrylovExt.jl +++ b/ext/UnfoldKrylovExt.jl @@ -1,75 +1,75 @@ module UnfoldKrylovExt - import Krylov - using CUDA, CUDA.CUSPARSE, SparseArrays - using Unfold - using Missings - using ProgressMeter +import Krylov +using CUDA, CUDA.CUSPARSE, SparseArrays +using Unfold +using Missings +using ProgressMeter - """ - Alternative implementation of LSMR using Krylov.jl +""" +Alternative implementation of LSMR using Krylov.jl - The fastest solver (x30 in some cases) with GPU = true +The fastest solver (x30 in some cases) with GPU = true - To activate you have to do: - > using Krylov,CUDA +To activate you have to do: +> using Krylov,CUDA - Difference to solver_default: No suport for per-channel missings. If one sample is missing in any channel, whole channel is removed due a lack of support for Missings in Krylov. +Difference to solver_default: No suport for per-channel missings. If one sample is missing in any channel, whole channel is removed due a lack of support for Missings in Krylov. - """ - function solver_krylov( - X, - data::AbstractArray{T,2}; - GPU=false, - history=true, - multithreading = GPU ? false : true, - showprogress=true, - stderror=false - ) where {T<:Union{Missing,<:Number}} +""" +function solver_krylov( + X, + data::AbstractArray{T,2}; + GPU = false, + history = true, + multithreading = GPU ? false : true, + showprogress = true, + stderror = false, +) where {T<:Union{Missing,<:Number}} @assert !(multithreading && GPU) "currently no support for both GPU and multi-threading" - minfo = Array{Any,1}(undef,size(data,1)) - - ix = any(@. !ismissing(data); dims=1)[1, :] - X_loop = disallowmissing(X[ix, :]) - data = disallowmissing(view(data, :, ix)) - - #@show typeof(X_loop) - if GPU - X_loop = CuSparseMatrixCSC(X_loop) - lsmr_solver = Krylov.LsmrSolver(size(X_loop)..., CuVector{Float64}) - data = CuArray(data) - else - lsmr_solver = Krylov.LsmrSolver(size(X_loop)..., Vector{Float64}) - end - - p = Progress(size(data,1);enabled=showprogress) - - beta = zeros(size(data, 1), size(X, 2)) # had issues with undef - Unfold.@maybe_threads multithreading for ch = 1:size(data, 1) - @debug ch - data_loop = view(data,ch,:) - # use the previous channel as a starting point - #ch == 1 || Krylov.warm_start!(lsmr_solver,beta[ch-1,:])#copyto!(view(beta, ch, :), view(beta, ch-1, :)) - - Krylov.solve!(lsmr_solver, X_loop, data_loop; history=history) - beta[ch, :] = Vector(lsmr_solver.x) - - #beta[ch,:],h = Krylov.lsmr(X_loop,data_loop,history=history) - - - minfo[ch] = deepcopy(lsmr_solver.stats) - next!(p) - end - finish!(p) - if stderror - stderror = Unfold.calculate_stderror(X, data, beta) - modelfit = Unfold.LinearModelFit(beta, ["krylov_lsmr", minfo], stderror) - else - modelfit = Unfold.LinearModelFit(beta, ["krylov_lsmr", minfo]) - end - return modelfit + minfo = Array{Any,1}(undef, size(data, 1)) + + ix = any(@. !ismissing(data); dims = 1)[1, :] + X_loop = disallowmissing(X[ix, :]) + data = disallowmissing(view(data, :, ix)) + + #@show typeof(X_loop) + if GPU + X_loop = CuSparseMatrixCSC(X_loop) + lsmr_solver = Krylov.LsmrSolver(size(X_loop)..., CuVector{Float64}) + data = CuArray(data) + else + lsmr_solver = Krylov.LsmrSolver(size(X_loop)..., Vector{Float64}) end + p = Progress(size(data, 1); enabled = showprogress) -end \ No newline at end of file + beta = zeros(size(data, 1), size(X, 2)) # had issues with undef + Unfold.@maybe_threads multithreading for ch = 1:size(data, 1) + @debug ch + data_loop = view(data, ch, :) + # use the previous channel as a starting point + #ch == 1 || Krylov.warm_start!(lsmr_solver,beta[ch-1,:])#copyto!(view(beta, ch, :), view(beta, ch-1, :)) + + Krylov.solve!(lsmr_solver, X_loop, data_loop; history = history) + beta[ch, :] = Vector(lsmr_solver.x) + + #beta[ch,:],h = Krylov.lsmr(X_loop,data_loop,history=history) + + + minfo[ch] = deepcopy(lsmr_solver.stats) + next!(p) + end + finish!(p) + if stderror + stderror = Unfold.calculate_stderror(X, data, beta) + modelfit = Unfold.LinearModelFit(beta, ["krylov_lsmr", minfo], stderror) + else + modelfit = Unfold.LinearModelFit(beta, ["krylov_lsmr", minfo]) + end + return modelfit +end + + +end diff --git a/ext/UnfoldMixedModelsExt/UnfoldMixedModelsExt.jl b/ext/UnfoldMixedModelsExt/UnfoldMixedModelsExt.jl index d0cff635..1d57b8d0 100644 --- a/ext/UnfoldMixedModelsExt/UnfoldMixedModelsExt.jl +++ b/ext/UnfoldMixedModelsExt/UnfoldMixedModelsExt.jl @@ -2,7 +2,7 @@ module UnfoldMixedModelsExt using Unfold -import Unfold: isMixedModelFormula,make_estimate +import Unfold: isMixedModelFormula, make_estimate using MixedModels #import MixedModels.FeMat # extended for sparse femats, type piracy => issue on MixedModels.jl github using StaticArrays # for MixedModels extraction of parametrs (inherited from MixedModels.jl, not strictly needed ) @@ -23,4 +23,4 @@ include("timeexpandedterm.jl") include("typedefinitions.jl") include("effects.jl") -end \ No newline at end of file +end diff --git a/ext/UnfoldMixedModelsExt/condense.jl b/ext/UnfoldMixedModelsExt/condense.jl index 9835a5cc..e69cdfe2 100644 --- a/ext/UnfoldMixedModelsExt/condense.jl +++ b/ext/UnfoldMixedModelsExt/condense.jl @@ -14,7 +14,7 @@ function Unfold.make_estimate( ) ranef_group = [x.group for x in MixedModels.tidyσs(modelfit(m))] - + # reshape to pred x time x chan and then invert to chan x time x pred ranef_group = permutedims( reshape(ranef_group, :, size(coef(m), 2), size(coef(m), 1)), @@ -35,7 +35,9 @@ function Unfold.make_estimate( return Float64.(estimate), stderror, group end -function Unfold.stderror(m::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime}) +function Unfold.stderror( + m::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime}, +) return permutedims( reshape(vcat([[b.se...] for b in modelfit(m).fits]...), reverse(size(coef(m)))), [3, 2, 1], @@ -45,13 +47,13 @@ end function Unfold.get_coefnames(uf::UnfoldLinearMixedModelContinuousTime) # special case here, because we have to reorder the random effects to the end, else labels get messed up as we concat (coefs,ranefs) - # coefnames = Unfold.coefnames(formula(uf)) -# coefnames(formula(uf)[1].rhs[1]) + # coefnames = Unfold.coefnames(formula(uf)) + # coefnames(formula(uf)[1].rhs[1]) formulas = Unfold.formula(uf) - if !isa(formulas,AbstractArray) # in case we have only a single basisfunction + if !isa(formulas, AbstractArray) # in case we have only a single basisfunction formulas = [formulas] end fe_coefnames = vcat([coefnames(f.rhs[1]) for f in formulas]...) re_coefnames = vcat([coefnames(f.rhs[2:end]) for f in formulas]...) - return vcat(fe_coefnames,re_coefnames) -end \ No newline at end of file + return vcat(fe_coefnames, re_coefnames) +end diff --git a/ext/UnfoldMixedModelsExt/designmatrix.jl b/ext/UnfoldMixedModelsExt/designmatrix.jl index 12cdcbc6..64ae1fbf 100644 --- a/ext/UnfoldMixedModelsExt/designmatrix.jl +++ b/ext/UnfoldMixedModelsExt/designmatrix.jl @@ -5,19 +5,23 @@ end check_groupsorting(r::MatrixTerm) = check_groupsorting(r.terms) function check_groupsorting(r::Tuple) -@debug "checking group sorting" -ix = findall(isa.(r,MixedModels.AbstractReTerm)) + @debug "checking group sorting" + ix = findall(isa.(r, MixedModels.AbstractReTerm)) -rhs(x::RandomEffectsTerm) = x.rhs -rhs(x::MixedModels.ZeroCorr) = rhs(x.term) -groupvars = [map(x->rhs(x).sym,r[ix])...] + rhs(x::RandomEffectsTerm) = x.rhs + rhs(x::MixedModels.ZeroCorr) = rhs(x.term) + groupvars = [map(x -> rhs(x).sym, r[ix])...] -@assert groupvars == sort(groupvars) "random effects have to be alphabetically ordered. e.g. (1+a|X) + (1+a|A) is not allowed. Please reorder" + @assert groupvars == sort(groupvars) "random effects have to be alphabetically ordered. e.g. (1+a|X) + (1+a|A) is not allowed. Please reorder" end -function Unfold.unfold_apply_schema(type::Type{<:Union{<:UnfoldLinearMixedModel,<:UnfoldLinearMixedModelContinuousTime}},f,schema) +function Unfold.unfold_apply_schema( + type::Type{<:Union{<:UnfoldLinearMixedModel,<:UnfoldLinearMixedModelContinuousTime}}, + f, + schema, +) @debug "LMM apply schema" - return apply_schema(f,schema, MixedModels.LinearMixedModel) + return apply_schema(f, schema, MixedModels.LinearMixedModel) end @@ -26,38 +30,38 @@ function StatsModels.coefnames(term::MixedModels.ZeroCorr) coefnames(term.term) end -function lmm_combineMats!(Xcomb,X1,X2) -# we have random effects - # combine REMats in single-eventtpe formulas ala y ~ (1|x) + (a|x) - Xs1 = MixedModels._amalgamate([X1.Xs[2:end]...], Float64) - Xs2 = MixedModels._amalgamate([X2.Xs[2:end]...], Float64) +function lmm_combineMats!(Xcomb, X1, X2) + # we have random effects + # combine REMats in single-eventtpe formulas ala y ~ (1|x) + (a|x) + Xs1 = MixedModels._amalgamate([X1.Xs[2:end]...], Float64) + Xs2 = MixedModels._amalgamate([X2.Xs[2:end]...], Float64) - Xcomb = (Xcomb, Xs1..., Xs2...) + Xcomb = (Xcomb, Xs1..., Xs2...) - # Next we make the ranefs all equal size - equalizeReMatLengths!(Xcomb[2:end]) + # Next we make the ranefs all equal size + equalizeReMatLengths!(Xcomb[2:end]) - # check if ranefs can be amalgamated. If this fails, then MixedModels tried to amalgamate over different eventtypes and we should throw the warning - # if it success, we have to check if the size before and after is identical. If it is not, it tried to amalgamize over different eventtypes which were of the same length + # check if ranefs can be amalgamated. If this fails, then MixedModels tried to amalgamate over different eventtypes and we should throw the warning + # if it success, we have to check if the size before and after is identical. If it is not, it tried to amalgamize over different eventtypes which were of the same length - try - reterms = MixedModels._amalgamate([Xcomb[2:end]...], Float64) - if length(reterms) != length(Xcomb[2:end]) - throw("IncompatibleRandomGroupings") - end - catch e - @error "Error, you seem to have two different eventtypes with the same random-effect grouping variable. \n - This is not allowed, you have to rename one. Example:\n - eventA: y~1+(1|item) \n - eventB: y~1+(1|item) \n - This leads to this error. Rename the later one\n - eventB: y~1+(1|itemB) " + try + reterms = MixedModels._amalgamate([Xcomb[2:end]...], Float64) + if length(reterms) != length(Xcomb[2:end]) throw("IncompatibleRandomGroupings") end + catch e + @error "Error, you seem to have two different eventtypes with the same random-effect grouping variable. \n + This is not allowed, you have to rename one. Example:\n + eventA: y~1+(1|item) \n + eventB: y~1+(1|item) \n + This leads to this error. Rename the later one\n + eventB: y~1+(1|itemB) " + throw("IncompatibleRandomGroupings") + end - return Xcomb + return Xcomb end - + function changeReMatSize!(remat::MixedModels.AbstractReMat, m::Integer) @@ -123,8 +127,10 @@ $(SIGNATURES) This function timeexpands the random effects and generates a ReMat object """ function StatsModels.modelcols( - term::Unfold.TimeExpandedTerm{<:Union{<:RandomEffectsTerm,<:MixedModels.AbstractReTerm}}, - tbl + term::Unfold.TimeExpandedTerm{ + <:Union{<:RandomEffectsTerm,<:MixedModels.AbstractReTerm}, + }, + tbl, ) # exchange this to get ZeroCorr to work @@ -132,10 +138,10 @@ function StatsModels.modelcols( tbl = DataFrame(tbl) # get the non-timeexpanded reMat reMat = modelcols(term.term, tbl) - + # Timeexpand the designmatrix z = transpose(Unfold.time_expand(transpose(reMat.z), term, tbl)) - + z = disallowmissing(z) # can't have missing in here @@ -245,6 +251,6 @@ function changeMatSize!(m, fe::AbstractSparseMatrix, remats) end function changeMatSize!(m, fe::AbstractMatrix, remats) changeReMatSize!.(remats, Ref(m)) - fe = fe[1:m,:] + fe = fe[1:m, :] return (fe, remats...) end diff --git a/ext/UnfoldMixedModelsExt/effects.jl b/ext/UnfoldMixedModelsExt/effects.jl index ef6d143b..2b054d39 100644 --- a/ext/UnfoldMixedModelsExt/effects.jl +++ b/ext/UnfoldMixedModelsExt/effects.jl @@ -1,4 +1,4 @@ function StatsModels.modelmatrix(model::UnfoldLinearMixedModel, bool) @assert bool == false "time continuous model matrix is not implemented for a `UnfoldLinearMixedModel`" return designmatrix(model).Xs[1] -end \ No newline at end of file +end diff --git a/ext/UnfoldMixedModelsExt/fit.jl b/ext/UnfoldMixedModelsExt/fit.jl index c8a8cafc..dafef10a 100644 --- a/ext/UnfoldMixedModelsExt/fit.jl +++ b/ext/UnfoldMixedModelsExt/fit.jl @@ -21,7 +21,8 @@ function StatsModels.fit!( #@assert length(first(values(design(uf)))[2]) if uf isa UnfoldLinearMixedModel if ~isempty(Unfold.design(uf)) - @assert length(Unfold.times(Unfold.design(uf))) == size(data,length(size(data))-1) "Times Vector does not match second last dimension of input data - forgot to epoch, or misspecified 'time' vector?" + @assert length(Unfold.times(Unfold.design(uf))) == + size(data, length(size(data)) - 1) "Times Vector does not match second last dimension of input data - forgot to epoch, or misspecified 'time' vector?" end end # function content partially taken from MixedModels.jl bootstrap.jl @@ -40,13 +41,13 @@ function StatsModels.fit!( end nchan = size(data, 1) - Xs = (Unfold.equalizeLengths(Xs[1]),Xs[2:end]...) - _,data = Unfold.zeropad(Xs[1],data) + Xs = (Unfold.equalizeLengths(Xs[1]), Xs[2:end]...) + _, data = Unfold.zeropad(Xs[1], data) # get a un-fitted mixed model object - + Xs = disallowmissing.(Xs) #Xs = (Matrix(Xs[1]),Xs[2:end]...) - + mm = LinearMixedModel_wrapper(Unfold.formula(uf), firstData, Xs) # prepare some variables to be used βsc, θsc = similar(MixedModels.coef(mm)), similar(mm.θ) # pre allocate @@ -134,7 +135,13 @@ function reshape_lmm(uf::UnfoldLinearMixedModelContinuousTime, est) end -LinearMixedModel_wrapper(form,data::Array{<:Union{TData},1},Xs;wts = []) where {TData<:Union{Missing,Number}}= @error("currently no support for missing values in MixedModels.jl") +LinearMixedModel_wrapper( + form, + data::Array{<:Union{TData},1}, + Xs; + wts = [], +) where {TData<:Union{Missing,Number}} = + @error("currently no support for missing values in MixedModels.jl") """ $(SIGNATURES) @@ -149,21 +156,21 @@ function LinearMixedModel_wrapper( wts = [], ) where {TData<:Number} # function LinearMixedModel_wrapper(form,data::Array{<:Union{Missing,TData},1},Xs;wts = []) where {TData<:Number} - Xs = (Unfold.equalizeLengths(Xs[1]),Xs[2:end]...) + Xs = (Unfold.equalizeLengths(Xs[1]), Xs[2:end]...) # XXX Push this to utilities zeropad # Make sure X & y are the same size m = size(Xs[1])[1] if m != size(data)[1] - fe,data = Unfold.zeropad(Xs[1],data) - + fe, data = Unfold.zeropad(Xs[1], data) + Xs = changeMatSize!(size(data)[1], fe, Xs[2:end]) end #y = (reshape(float(data), (:, 1))) y = data - + MixedModels.LinearMixedModel(y, Xs, form, wts) end @@ -183,4 +190,4 @@ function MixedModels.LinearMixedModel(y, Xs, form::Array, wts) end -isMixedModelFormula(f::typeof(MixedModels.zerocorr)) = true \ No newline at end of file +isMixedModelFormula(f::typeof(MixedModels.zerocorr)) = true diff --git a/ext/UnfoldMixedModelsExt/statistics.jl b/ext/UnfoldMixedModelsExt/statistics.jl index 5b7b8269..4b68ab90 100644 --- a/ext/UnfoldMixedModelsExt/statistics.jl +++ b/ext/UnfoldMixedModelsExt/statistics.jl @@ -4,10 +4,10 @@ $(SIGNATURES) Returns a partial LMM model (non-functional due to lacking data) to be used in likelihoodratiotests. `k` to selcet which of the modelfit's to fake """ -function fake_lmm(m::UnfoldLinearMixedModel,k::Int) - (feterm,reterm) = Unfold.designmatrix(m).Xs - fakeY = zeros(size(feterm,1)) - lmm = LinearMixedModel_wrapper(Unfold.formula(m),fakeY,designmatrix(m).Xs) +function fake_lmm(m::UnfoldLinearMixedModel, k::Int) + (feterm, reterm) = Unfold.designmatrix(m).Xs + fakeY = zeros(size(feterm, 1)) + lmm = LinearMixedModel_wrapper(Unfold.formula(m), fakeY, designmatrix(m).Xs) fcoll = Unfold.modelfit(m) #lmm.objective .= fcoll.fits[1].objective @@ -16,62 +16,67 @@ function fake_lmm(m::UnfoldLinearMixedModel,k::Int) lmm.optsum.sigma = fcoll.fits[k].σ lmm.optsum.optimizer = :unfold lmm.optsum.returnvalue = :success - setθ!(lmm,fcoll.fits[k].θ) + setθ!(lmm, fcoll.fits[k].θ) return lmm end - -fake_lmm(m::UnfoldLinearMixedModel) = fake_lmm.(m,1:length(modelfit(m).fits)) - - """ - $(SIGNATURES) - - Calculate likelihoodratiotest - """ - function MixedModels.likelihoodratiotest(m::UnfoldLinearMixedModel...) - #@info lrtest(fake_lmm.(m,1)...) - n = length(Unfold.modelfit(m[1]).fits) - lrt = Array{MixedModels.LikelihoodRatioTest}(undef,n) - sizehint!(lrt,n) - for k = 1:n - ms = fake_lmm.(m,k) - #@info objective.(ms) - lrt[k] = MixedModels.likelihoodratiotest(ms...) - end - return lrt - end - - """ - $(SIGNATURES) - Unfold-Method: return pvalues of likelihoodratiotests, typically calculated: - - # Examples - julia> pvalues(likelihoodratiotest(m1,m2)) - - where m1/m2 are UnfoldLinearMixedModel's - - Tipp: if you only compare two models you can easily get a vector of p-values: - - julia> vcat(pvalues(likelihoodratiotest(m1,m2))...) - - - Multiple channels are returned linearized at the moment, as we do not have access to the amount of channels after the LRT, you can do: - - julia> reshape(vcat(pvalues(likelihoodratiotest(m1,m2))...),ntimes,nchan)' - - """ - function pvalues(lrtvec::Vector{MixedModels.LikelihoodRatioTest}) - [lrt.pvalues for lrt in lrtvec] - end - - - - function MixedModels.rePCA(m::UnfoldLinearMixedModel) - oneresult = MixedModels.rePCA(fake_lmm(m,1)) - emptyPCA = (;) - for s = keys(oneresult) - res = hcat([a[s] for a in MixedModels.rePCA.(fake_lmm.(Ref(m),1:length(modelfit(m).fits)))]...) - newPCA = NamedTuple{(s,)}([res]) - emptyPCA = merge(emptyPCA,newPCA) - end - return emptyPCA - end \ No newline at end of file + +fake_lmm(m::UnfoldLinearMixedModel) = fake_lmm.(m, 1:length(modelfit(m).fits)) + +""" +$(SIGNATURES) + +Calculate likelihoodratiotest +""" +function MixedModels.likelihoodratiotest(m::UnfoldLinearMixedModel...) + #@info lrtest(fake_lmm.(m,1)...) + n = length(Unfold.modelfit(m[1]).fits) + lrt = Array{MixedModels.LikelihoodRatioTest}(undef, n) + sizehint!(lrt, n) + for k = 1:n + ms = fake_lmm.(m, k) + #@info objective.(ms) + lrt[k] = MixedModels.likelihoodratiotest(ms...) + end + return lrt +end + +""" +$(SIGNATURES) +Unfold-Method: return pvalues of likelihoodratiotests, typically calculated: + +# Examples +julia> pvalues(likelihoodratiotest(m1,m2)) + +where m1/m2 are UnfoldLinearMixedModel's + +Tipp: if you only compare two models you can easily get a vector of p-values: + +julia> vcat(pvalues(likelihoodratiotest(m1,m2))...) + + +Multiple channels are returned linearized at the moment, as we do not have access to the amount of channels after the LRT, you can do: + +julia> reshape(vcat(pvalues(likelihoodratiotest(m1,m2))...),ntimes,nchan)' + +""" +function pvalues(lrtvec::Vector{MixedModels.LikelihoodRatioTest}) + [lrt.pvalues for lrt in lrtvec] +end + + + +function MixedModels.rePCA(m::UnfoldLinearMixedModel) + oneresult = MixedModels.rePCA(fake_lmm(m, 1)) + emptyPCA = (;) + for s in keys(oneresult) + res = hcat( + [ + a[s] for + a in MixedModels.rePCA.(fake_lmm.(Ref(m), 1:length(modelfit(m).fits))) + ]..., + ) + newPCA = NamedTuple{(s,)}([res]) + emptyPCA = merge(emptyPCA, newPCA) + end + return emptyPCA +end diff --git a/ext/UnfoldMixedModelsExt/timeexpandedterm.jl b/ext/UnfoldMixedModelsExt/timeexpandedterm.jl index c29f21c1..9ed8f377 100644 --- a/ext/UnfoldMixedModelsExt/timeexpandedterm.jl +++ b/ext/UnfoldMixedModelsExt/timeexpandedterm.jl @@ -6,4 +6,4 @@ function Unfold.TimeExpandedTerm( ) where {N} # for mixed models, apply it to each Term Unfold.TimeExpandedTerm.(term, Ref(basisfunction), Ref(eventfields)) -end \ No newline at end of file +end diff --git a/ext/UnfoldMixedModelsExt/typedefinitions.jl b/ext/UnfoldMixedModelsExt/typedefinitions.jl index 132bcb3c..b861410e 100644 --- a/ext/UnfoldMixedModelsExt/typedefinitions.jl +++ b/ext/UnfoldMixedModelsExt/typedefinitions.jl @@ -1,9 +1,9 @@ struct UnfoldMixedModelFitCollection{T<:AbstractFloat} <: - MixedModels.MixedModelFitCollection{T} - fits::Vector - λ::Vector{<:Union{LowerTriangular{T,Matrix{T}},Diagonal{T,Vector{T}}}} - inds::Vector{Vector{Int}} - lowerbd::Vector{T} - fcnames::NamedTuple -end \ No newline at end of file + MixedModels.MixedModelFitCollection{T} + fits::Vector + λ::Vector{<:Union{LowerTriangular{T,Matrix{T}},Diagonal{T,Vector{T}}}} + inds::Vector{Vector{Int}} + lowerbd::Vector{T} + fcnames::NamedTuple +end diff --git a/ext/UnfoldRobustModelsExt.jl b/ext/UnfoldRobustModelsExt.jl index 3a492a11..9b57cbb3 100644 --- a/ext/UnfoldRobustModelsExt.jl +++ b/ext/UnfoldRobustModelsExt.jl @@ -8,44 +8,44 @@ using Missings function solver_robust( X, data::AbstractArray{T,3}; - estimator = MEstimator{TukeyLoss}(), - rlmOptions = (initial_scale=:mad,) + estimator = MEstimator{TukeyLoss}(), + rlmOptions = (initial_scale = :mad,), ) where {T<:Union{Missing,<:Number}} #beta = zeros(Union{Missing,Number},size(data, 1), size(data, 2), size(X, 2)) - beta = zeros(T,size(data, 1), size(data, 2), size(X, 2)) + beta = zeros(T, size(data, 1), size(data, 2), size(X, 2)) @showprogress 0.1 for ch = 1:size(data, 1) for t = 1:size(data, 2) #@debug("$(ndims(data,)),$t,$ch") - - dd = view(data, ch, t,:) + + dd = view(data, ch, t, :) ix = @. !ismissing(dd) - # init beta - #if ch == 1 && t==1 - #elseif ch > 1 && t == 1 - # copyto!(view(beta, ch, 1,:), view(beta, ch-1, 1,:)) - #else - # copyto!(view(beta, ch,t, :), view(beta, ch,t-1, :)) - #end - - X_local = disallowmissing((X[ix,:])) # view crashes robust model here. XXX follow up once + # init beta + #if ch == 1 && t==1 + #elseif ch > 1 && t == 1 + # copyto!(view(beta, ch, 1,:), view(beta, ch-1, 1,:)) + #else + # copyto!(view(beta, ch,t, :), view(beta, ch,t-1, :)) + #end + + X_local = disallowmissing((X[ix, :])) # view crashes robust model here. XXX follow up once # https://github.com/JuliaStats/GLM.jl/issues/470 received a satisfying result - y_local = disallowmissing(@view(data[ch,t,ix])) + y_local = disallowmissing(@view(data[ch, t, ix])) # if not specified otherwise # this results in strange regularizing effects - #if :initial_coef ∉ keys(rlmOptions) + #if :initial_coef ∉ keys(rlmOptions) # rlmOptions = merge(rlmOptions,(initial_coef=@view(beta[ch,t,:]),)) #end - m = rlm(X_local,y_local,estimator;rlmOptions...) - beta[ch,t,:] .= coef(m) - + m = rlm(X_local, y_local, estimator; rlmOptions...) + beta[ch, t, :] .= coef(m) + end end - modelfit = Unfold.LinearModelFit(beta, ["solver_robust"]) - + modelfit = Unfold.LinearModelFit(beta, ["solver_robust"]) + return modelfit end -end \ No newline at end of file +end diff --git a/src/Unfold.jl b/src/Unfold.jl index bf66b694..1df7371c 100644 --- a/src/Unfold.jl +++ b/src/Unfold.jl @@ -96,73 +96,73 @@ if !isdefined(Base, :get_extension) solver_krylov = UnfoldKrylovExt.solver_krylov else # Jl 1.9+: we need dummy functions, in case the extension was not activated to warn the user if they try to use a functionality that is not yet defined - checkFun(sym) = Base.get_extension(@__MODULE__(),sym) - function solver_robust(args...;kwargs...) - ext = checkFun(:UnfoldRobustModelsExt) + checkFun(sym) = Base.get_extension(@__MODULE__(), sym) + function solver_robust(args...; kwargs...) + ext = checkFun(:UnfoldRobustModelsExt) msg = "RobustModels not loaded. Please use ]add RobustModels, using RobustModels to install it prior to using" - isnothing(ext) ? throw(msg) : ext.solver_robust(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.solver_robust(args...; kwargs...) end - function pvalues(args...;kwargs...) - ext = checkFun(:UnfoldMixedModelsExt) + function pvalues(args...; kwargs...) + ext = checkFun(:UnfoldMixedModelsExt) msg = "MixedModels not loaded. Please use ]add MixedModels, using MixedModels to install it prior to using" - isnothing(ext) ? throw(msg) : ext.pvalues(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.pvalues(args...; kwargs...) end - function likelihoodratiotest(args...;kwargs...) - ext = checkFun(:UnfoldMixedModelsExt) + function likelihoodratiotest(args...; kwargs...) + ext = checkFun(:UnfoldMixedModelsExt) msg = "MixedModels not loaded. Please use ]add MixedModels, using MixedModels to install it prior to using" - isnothing(ext) ? throw(msg) : ext.likelihoodratiotest(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.likelihoodratiotest(args...; kwargs...) end - function rePCA(args...;kwargs...) - ext = checkFun(:UnfoldMixedModelsExt) + function rePCA(args...; kwargs...) + ext = checkFun(:UnfoldMixedModelsExt) msg = "MixedModels not loaded. Please use ]add MixedModels, using MixedModels to install it prior to using" - isnothing(ext) ? throw(msg) : ext.rePCA(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.rePCA(args...; kwargs...) end - function check_groupsorting(args...;kwargs...) - ext = checkFun(:UnfoldMixedModelsExt) + function check_groupsorting(args...; kwargs...) + ext = checkFun(:UnfoldMixedModelsExt) msg = "MixedModels not loaded. Please use ]add MixedModels, using MixedModels to install it prior to using" - isnothing(ext) ? throw(msg) : ext.check_groupsorting(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.check_groupsorting(args...; kwargs...) end - function lmm_combineMats!(args...;kwargs...) - ext = checkFun(:UnfoldMixedModelsExt) + function lmm_combineMats!(args...; kwargs...) + ext = checkFun(:UnfoldMixedModelsExt) msg = "MixedModels not loaded. Please use ]add MixedModels, using MixedModels to install it prior to using" - isnothing(ext) ? throw(msg) : ext.lmm_combineMats!(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.lmm_combineMats!(args...; kwargs...) end - function splinebasis(args...;kwargs...) - ext = checkFun(:UnfoldBSplineKitExt) + function splinebasis(args...; kwargs...) + ext = checkFun(:UnfoldBSplineKitExt) msg = "BSplineKit not loaded. Please use `]add BSplineKit, using BSplineKit` to install/load it, if you want to use splines" - isnothing(ext) ? throw(msg) : ext.splinebasis(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.splinebasis(args...; kwargs...) end - - - - function spl(args...;kwargs...) - ext = checkFun(:UnfoldBSplineKitExt) + + + + function spl(args...; kwargs...) + ext = checkFun(:UnfoldBSplineKitExt) msg = "BSplineKit not loaded. Please use `]add BSplineKit, using BSplineKit` to install/load it, if you want to use splines" - isnothing(ext) ? throw(msg) : ext.spl(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.spl(args...; kwargs...) end - function circspl(args...;kwargs...) - ext = checkFun(:UnfoldBSplineKitExt) + function circspl(args...; kwargs...) + ext = checkFun(:UnfoldBSplineKitExt) msg = "BSplineKit not loaded. Please use `]add BSplineKit, using BSplineKit` to install/load it, if you want to use splines" - isnothing(ext) ? throw(msg) : ext.circspl(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.circspl(args...; kwargs...) end - - function solver_krylov(args...;kwargs...) - ext = checkFun(:UnfoldKrylovExt) + + function solver_krylov(args...; kwargs...) + ext = checkFun(:UnfoldKrylovExt) msg = "Krylov or CUDA not loaded. Please use `]add Krylov,CUDA, using Krylov,CUDA` to install/load it, if you want to use GPU-fitting" - isnothing(ext) ? throw(msg) : ext.solver_krylov(args...;kwargs...) + isnothing(ext) ? throw(msg) : ext.solver_krylov(args...; kwargs...) end - + end - -export spl,circspl + +export spl, circspl export likelihoodratiotest # statistics.jl export pvalues # statistics.jl diff --git a/src/basisfunctions.jl b/src/basisfunctions.jl index 731f59b0..b88b5652 100644 --- a/src/basisfunctions.jl +++ b/src/basisfunctions.jl @@ -37,7 +37,11 @@ struct FIRBasis <: BasisFunction shiftOnset::Int64 end -@deprecate FIRBasis(kernel::Function,times,name,shiftOnset) FIRBasis(times,name,shiftOnset) +@deprecate FIRBasis(kernel::Function, times, name, shiftOnset) FIRBasis( + times, + name, + shiftOnset, +) collabel(basis::FIRBasis) = :time colnames(basis::FIRBasis) = basis.times[1:end-1] @@ -98,7 +102,7 @@ julia> f(103.3) """ function firbasis(τ, sfreq, name::String) τ = round_times(τ, sfreq) - times = range(τ[1], stop = τ[2]+1 ./ sfreq, step = 1 ./ sfreq) # stop + 1 step, because we support fractional event-timings + times = range(τ[1], stop = τ[2] + 1 ./ sfreq, step = 1 ./ sfreq) # stop + 1 step, because we support fractional event-timings @@ -131,9 +135,13 @@ function firkernel(e, times) eboth[isapprox.(eboth, 0, atol = 1e-15)] .= 0 ksize = length(times) # kernelsize - kernel = - spdiagm(ksize + 1, ksize, 0 => repeat([eboth[1]], ksize), -1 => repeat([eboth[2]], ksize)) - + kernel = spdiagm( + ksize + 1, + ksize, + 0 => repeat([eboth[1]], ksize), + -1 => repeat([eboth[2]], ksize), + ) + # dropzeros!(kernel) # we often get implicit 0, especially if the latencies are rounded return (kernel) @@ -185,7 +193,12 @@ function hrfbasis( # p(6) - onset {seconds} 0 # p(7) - length of kernel {seconds} 32 kernel = e -> hrfkernel(e, TR, parameters) - return HRFBasis(kernel,["f(x)"],range(0, (length(kernel([0, 1])) - 1) * TR, step = TR),name) + return HRFBasis( + kernel, + ["f(x)"], + range(0, (length(kernel([0, 1])) - 1) * TR, step = TR), + name, + ) end shiftOnset(basis::HRFBasis) = 0 @@ -203,10 +216,10 @@ collabel(term::Array{<:AbstractTerm}) = collabel(term[1].rhs) # in case of comb shiftOnset(basis::BasisFunction) = basis.shiftOnset colnames(basis::BasisFunction) = basis.colnames -kernel(basis::BasisFunction,e) = basis.kernel(e) +kernel(basis::BasisFunction, e) = basis.kernel(e) @deprecate kernel(basis::BasisFunction) basis.kernel -kernel(basis::FIRBasis,e) = firkernel(e,basis.times[1:end-1]) +kernel(basis::FIRBasis, e) = firkernel(e, basis.times[1:end-1]) times(basis::BasisFunction) = basis.times name(basis::BasisFunction) = basis.name @@ -268,29 +281,29 @@ end # taken from https://codereview.stackexchange.com/questions/284537/implementing-a-1d-convolution-simd-friendly-in-julia # to replace the DSP.conv function -function conv1D!( vO :: Array{T, 1}, vA :: Array{T, 1}, vB :: Array{T, 1} ) :: Array{T, 1} where {T <: Real} +function conv1D!(vO::Array{T,1}, vA::Array{T,1}, vB::Array{T,1})::Array{T,1} where {T<:Real} - lenA = length(vA); - lenB = length(vB); + lenA = length(vA) + lenB = length(vB) - fill!(vO, zero(T)); - for idxB in 1:lenB - for idxA in 1:lenA - @inbounds vO[idxA + idxB - 1] += vA[idxA] * vB[idxB]; + fill!(vO, zero(T)) + for idxB = 1:lenB + for idxA = 1:lenA + @inbounds vO[idxA+idxB-1] += vA[idxA] * vB[idxB] end end - return vO; + return vO end -function conv1D( vA :: Array{T, 1}, vB :: Array{T, 1} ) :: Array{T, 1} where {T <: Real} +function conv1D(vA::Array{T,1}, vB::Array{T,1})::Array{T,1} where {T<:Real} - lenA = length(vA); - lenB = length(vB); + lenA = length(vA) + lenB = length(vB) - vO = Array{T, 1}(undef, lenA + lenB - 1); + vO = Array{T,1}(undef, lenA + lenB - 1) - return conv1D!(vO, vA, vB); + return conv1D!(vO, vA, vB) -end \ No newline at end of file +end diff --git a/src/condense.jl b/src/condense.jl index 6471934f..d88d9d53 100644 --- a/src/condense.jl +++ b/src/condense.jl @@ -11,7 +11,7 @@ function get_coefnames(uf::Union{UnfoldModel,DesignMatrix}) end - + modelfit(uf::UnfoldModel) = uf.modelfit StatsModels.coef(uf::UnfoldModel) = coef(modelfit(uf)) StatsModels.coef(mf::LinearModelFit) = mf.estimate @@ -80,22 +80,22 @@ function StatsModels.coeftable(uf::Union{UnfoldLinearModel,UnfoldLinearMixedMode colnames_basis_rep = permutedims(repeat(colnames_basis, 1, nchan, ncoefs), [2 1 3]) chan_rep = repeat(1:nchan, 1, ncols, ncoefs) - designkeys = collect(keys(design(uf))) + designkeys = collect(keys(design(uf))) + + + - - - if length(designkeys) == 1 # in case of 1 event, repeat it by ncoefs - basisnames = repeat(["event: $(designkeys[1])"],ncoefs) + basisnames = repeat(["event: $(designkeys[1])"], ncoefs) else basisnames = String[] - for (ix,evt) = enumerate(designkeys) - push!(basisnames,repeat(["event: $(evt)"],size(modelmatrix(uf)[ix],2) )...) + for (ix, evt) in enumerate(designkeys) + push!(basisnames, repeat(["event: $(evt)"], size(modelmatrix(uf)[ix], 2))...) end end - - basisnames_rep = permutedims(repeat(basisnames,1, nchan, ncols),[2,3,1]) + + basisnames_rep = permutedims(repeat(basisnames, 1, nchan, ncols), [2, 3, 1]) # results = make_long_df(uf, coefs_rep, chan_rep, colnames_basis_rep, basisnames_rep, :time) diff --git a/src/designmatrix.jl b/src/designmatrix.jl index 76cf4e88..ce36236b 100644 --- a/src/designmatrix.jl +++ b/src/designmatrix.jl @@ -19,28 +19,28 @@ julia> Xdc = Xdc1+Xdc2 function combineDesignmatrices(X1::DesignMatrix, X2::DesignMatrix) # the reason for the assertion is simply that I found it too difficult to concatenate the formulas down below ;) should be easy to implement hough - @assert !(isa(X1.formulas,AbstractArray) && isa(X2.formulas,AbstractArray)) "it is currently not possible to combine desigmatrices from two already concatenated designs - please concatenate one after the other" + @assert !(isa(X1.formulas, AbstractArray) && isa(X2.formulas, AbstractArray)) "it is currently not possible to combine desigmatrices from two already concatenated designs - please concatenate one after the other" X1 = deepcopy(X1) X2 = deepcopy(X2) Xs1 = get_Xs(X1) Xs2 = get_Xs(X2) - if typeof(Xs1) <:Vector - Xcomb = [Xs1...,Xs2] + if typeof(Xs1) <: Vector + Xcomb = [Xs1..., Xs2] else - Xcomb = [Xs1,Xs2] + Xcomb = [Xs1, Xs2] end if typeof(X1.Xs) <: Tuple - Xcomb = lmm_combineMats!(Xcomb,X1,X2) + Xcomb = lmm_combineMats!(Xcomb, X1, X2) end if X1.formulas isa FormulaTerm # due to the assertion above, we can assume we have only 2 formulas here if X1.formulas.rhs isa Unfold.TimeExpandedTerm - fcomb = Vector{FormulaTerm{<:InterceptTerm, <:TimeExpandedTerm}}(undef,2) + fcomb = Vector{FormulaTerm{<:InterceptTerm,<:TimeExpandedTerm}}(undef, 2) else - fcomb = Vector{FormulaTerm}(undef,2) # mass univariate case + fcomb = Vector{FormulaTerm}(undef, 2) # mass univariate case end fcomb[1] = X1.formulas fcomb[2] = X2.formulas @@ -48,9 +48,12 @@ function combineDesignmatrices(X1::DesignMatrix, X2::DesignMatrix) else if X1.formulas[1].rhs isa Unfold.TimeExpandedTerm # we can ignore length of X2, as it has to be a single formula due to the assertion above - fcomb = Vector{FormulaTerm{<:InterceptTerm, <:TimeExpandedTerm}}(undef,length(X1.formulas)+1) + fcomb = Vector{FormulaTerm{<:InterceptTerm,<:TimeExpandedTerm}}( + undef, + length(X1.formulas) + 1, + ) else - fcomb = Vector{FormulaTerm}(undef,length(X1.formulas)+1) # mass univariate case + fcomb = Vector{FormulaTerm}(undef, length(X1.formulas) + 1) # mass univariate case end fcomb[1:end-1] = X1.formulas fcomb[end] = X2.formulas @@ -115,28 +118,31 @@ function designmatrix( @debug("generating DesignMatrix") # check for missings-columns - currently missings not really supported in StatsModels - s_tmp = schema(f,tbl,contrasts) # temporary scheme to get necessary terms + s_tmp = schema(f, tbl, contrasts) # temporary scheme to get necessary terms neededCols = [v.sym for v in values(s_tmp.schema)] tbl_nomissing = DataFrame(tbl) # tbl might be a SubDataFrame due to a view - but we can't modify a subdataframe, can we? - try - disallowmissing!(tbl_nomissing,neededCols) # if this fails, it means we have a missing value in a column we need. We do not support this + try + disallowmissing!(tbl_nomissing, neededCols) # if this fails, it means we have a missing value in a column we need. We do not support this catch e if e isa ArgumentError - error(e.msg*"\n we tried to get rid of a event-column declared as type Union{Missing,T}. But there seems to be some actual missing values in there. - You have to replace them yourself (e.g. replace(tbl.colWithMissingValue,missing=>0)) or impute them otherwise.") + error( + e.msg * + "\n we tried to get rid of a event-column declared as type Union{Missing,T}. But there seems to be some actual missing values in there. + You have to replace them yourself (e.g. replace(tbl.colWithMissingValue,missing=>0)) or impute them otherwise.", + ) else - rethrow() + rethrow() end end - - form = unfold_apply_schema(type,f, schema(f, tbl_nomissing, contrasts)) - + + form = unfold_apply_schema(type, f, schema(f, tbl_nomissing, contrasts)) + @debug "type: $type" - if (type== UnfoldLinearMixedModel) || (type == UnfoldLinearMixedModelContinuousTime) + if (type == UnfoldLinearMixedModel) || (type == UnfoldLinearMixedModelContinuousTime) # get all random effects check_groupsorting(form.rhs) @@ -145,7 +151,7 @@ function designmatrix( apply_basisfunction(form, basisfunction, get(Dict(kwargs), :eventfields, nothing)) # Evaluate the designmatrix - + #note that we use tbl again, not tbl_nomissing. @debug typeof(form) if (!isnothing(basisfunction)) & (type <: UnfoldLinearMixedModel) @@ -165,7 +171,7 @@ wrapper to make apply_schema mixed models as extension possible Note: type is not necessary here, but for LMM it is for multiple dispatch reasons! """ -unfold_apply_schema(type,f,schema) = apply_schema(f,schema,UnfoldModel) +unfold_apply_schema(type, f, schema) = apply_schema(f, schema, UnfoldModel) # specify for abstract interface @@ -181,8 +187,8 @@ function designmatrix( X = nothing fDict = design(uf) for (eventname, f) in pairs(fDict) - - @debug "Eventname, X:",eventname,X + + @debug "Eventname, X:", eventname, X if eventname == Any eventTbl = tbl else @@ -261,7 +267,7 @@ function designmatrix!(uf::UnfoldModel, evts; kwargs...) return uf end -function StatsModels.modelmatrix(uf::UnfoldLinearModel,basisfunction) +function StatsModels.modelmatrix(uf::UnfoldLinearModel, basisfunction) if basisfunction @warn("basisfunction not defined for this kind of model") else @@ -273,15 +279,16 @@ end equalizeLengths(Xs::AbstractMatrix) = Xs # UnfoldLinearMixedModelContinuousTime case -equalizeLengths(Xs::Tuple) = (equalizeLengths(Xs[1]),Xs[2:end]...) +equalizeLengths(Xs::Tuple) = (equalizeLengths(Xs[1]), Xs[2:end]...) # UnfoldLinearModel - they have to be equal already -equalizeLengths(Xs::Vector{<:AbstractMatrix}) = Xs +equalizeLengths(Xs::Vector{<:AbstractMatrix}) = Xs #UnfoldLinearModelContinuousTime -equalizeLengths(Xs::Vector{<:SparseMatrixCSC}) = equalizeLengths(Xs...) -equalizeLengths(Xs1::SparseMatrixCSC,Xs2::SparseMatrixCSC,args...) = equalizeLengths(equalizeLengths(Xs1,Xs2),args...) -function equalizeLengths(Xs1::SparseMatrixCSC,Xs2::SparseMatrixCSC) +equalizeLengths(Xs::Vector{<:SparseMatrixCSC}) = equalizeLengths(Xs...) +equalizeLengths(Xs1::SparseMatrixCSC, Xs2::SparseMatrixCSC, args...) = + equalizeLengths(equalizeLengths(Xs1, Xs2), args...) +function equalizeLengths(Xs1::SparseMatrixCSC, Xs2::SparseMatrixCSC) sX1 = size(Xs1, 1) sX2 = size(Xs2, 1) @@ -291,27 +298,27 @@ function equalizeLengths(Xs1::SparseMatrixCSC,Xs2::SparseMatrixCSC) elseif sX2 < sX1 Xs2 = SparseMatrixCSC(sX1, Xs2.n, Xs2.colptr, Xs2.rowval, Xs2.nzval) end - return hcat(Xs1,Xs2) + return hcat(Xs1, Xs2) end -function StatsModels.modelmatrix(uf::UnfoldLinearModelContinuousTime,basisfunction = true) +function StatsModels.modelmatrix(uf::UnfoldLinearModelContinuousTime, basisfunction = true) if basisfunction return modelmatrix(designmatrix(uf)) #return hcat(modelmatrix(designmatrix(uf))...) else # replace basisfunction with non-timeexpanded one f = formula(uf) - + # probably a more julian way to do this... - if isa(f,AbstractArray) - return modelcols_nobasis.(f,events(uf)) + if isa(f, AbstractArray) + return modelcols_nobasis.(f, events(uf)) else - return modelcols_nobasis(f,events(uf)) + return modelcols_nobasis(f, events(uf)) end - + end end -modelcols_nobasis(f::FormulaTerm,tbl::AbstractDataFrame) = modelcols(f.rhs.term,tbl) +modelcols_nobasis(f::FormulaTerm, tbl::AbstractDataFrame) = modelcols(f.rhs.term, tbl) StatsModels.modelmatrix(uf::UnfoldModel) = modelmatrix(designmatrix(uf))#modelmatrix(uf.design,uf.designmatrix.events) StatsModels.modelmatrix(d::DesignMatrix) = equalizeLengths(d.Xs) @@ -333,7 +340,7 @@ function formula(d::Dict) #TODO Specify Dict better else return [c[1] for c in collect(values(d))] end - + end """ @@ -342,28 +349,28 @@ calculates in which rows the individual event-basisfunctions should go in Xdc timeexpand_rows timeexpand_vals """ -function timeexpand_rows(onsets,bases,shift,ncolsX) +function timeexpand_rows(onsets, bases, shift, ncolsX) # generate rowindices rows = copy(rowvals.(bases)) - + # this shift is necessary as some basisfunction time-points can be negative. But a matrix is always from 1:τ. Thus we have to shift it backwards in time. # The onsets are onsets-1 XXX not sure why. - for r = eachindex(rows) + for r in eachindex(rows) rows[r] .+= floor(onsets[r] - 1) .+ shift end - - rows_red = reduce(vcat,rows) + + rows_red = reduce(vcat, rows) rows_red = repeat(rows_red, ncolsX) return rows_red - end +end """ $(SIGNATURES) calculates the actual designmatrix for a timeexpandedterm. Multiple dispatch on StatsModels.modelcols """ function StatsModels.modelcols(term::TimeExpandedTerm, tbl) - @debug term.term , first(tbl) + @debug term.term, first(tbl) X = modelcols(term.term, tbl) time_expand(X, term, tbl) @@ -384,13 +391,15 @@ function time_expand_getTimeRange(onset, basisfunction) end -function timeexpand_cols_allsamecols(bases,ncolsBasis::Int,ncolsX) +function timeexpand_cols_allsamecols(bases, ncolsBasis::Int, ncolsX) repeatEach = length(nzrange(bases[1], 1)) - cols_r = UnitRange{Int64}[((1:ncolsBasis) .+ix*ncolsBasis) for ix in (0:ncolsX-1) for b in 1:length(bases)] - - cols = reduce(vcat,cols_r) - - cols = repeat(cols,inner=repeatEach) + cols_r = UnitRange{Int64}[ + ((1:ncolsBasis) .+ ix * ncolsBasis) for ix in (0:ncolsX-1) for b = 1:length(bases) + ] + + cols = reduce(vcat, cols_r) + + cols = repeat(cols, inner = repeatEach) return cols end @@ -402,55 +411,53 @@ calculates in which rows the individual event-basisfunctions should go in Xdc see also timeexpand_rows timeexpand_vals """ -function timeexpand_cols(term,bases,ncolsBasis,ncolsX) +function timeexpand_cols(term, bases, ncolsBasis, ncolsX) # we can generate the columns much faster, if all bases output the same number of columns - fastpath = time_expand_allBasesSameCols(term.basisfunction,bases,ncolsBasis) + fastpath = time_expand_allBasesSameCols(term.basisfunction, bases, ncolsBasis) if fastpath - return timeexpand_cols_allsamecols(bases,ncolsBasis,ncolsX) + return timeexpand_cols_allsamecols(bases, ncolsBasis, ncolsX) else - return timeexpand_cols_generic(bases,ncolsBasis,ncolsX) + return timeexpand_cols_generic(bases, ncolsBasis, ncolsX) end end -function timeexpand_cols_generic(bases,ncolsBasis,ncolsX) - # it could happen, e.g. for bases that are duration modulated, that each event has different amount of columns - # in that case, we have to go the slow route - cols = Vector{Int64}[] - - for Xcol = 1:ncolsX - for b = 1:length(bases) - coloffset = (Xcol - 1) * ncolsBasis - for c = 1:ncolsBasis - push!(cols, - repeat([c + coloffset], length(nzrange(bases[b], c))), - ) - end +function timeexpand_cols_generic(bases, ncolsBasis, ncolsX) + # it could happen, e.g. for bases that are duration modulated, that each event has different amount of columns + # in that case, we have to go the slow route + cols = Vector{Int64}[] + + for Xcol = 1:ncolsX + for b = 1:length(bases) + coloffset = (Xcol - 1) * ncolsBasis + for c = 1:ncolsBasis + push!(cols, repeat([c + coloffset], length(nzrange(bases[b], c)))) end end - return reduce(vcat,cols) - + end + return reduce(vcat, cols) + end -function timeexpand_vals(bases,X,nTotal,ncolsX) - # generate values - #vals = [] - vals = Array{Union{Missing,Float64}}(undef,nTotal) - ix = 1 - - for Xcol = 1:ncolsX - for (i,b) = enumerate(bases) - b_nz = nonzeros(b) - l = length(b_nz) - - vals[ix:ix+l-1] .= b_nz .* @view X[i, Xcol] - ix = ix+l +function timeexpand_vals(bases, X, nTotal, ncolsX) + # generate values + #vals = [] + vals = Array{Union{Missing,Float64}}(undef, nTotal) + ix = 1 + + for Xcol = 1:ncolsX + for (i, b) in enumerate(bases) + b_nz = nonzeros(b) + l = length(b_nz) + + vals[ix:ix+l-1] .= b_nz .* @view X[i, Xcol] + ix = ix + l #push!(vals, ) - end end + end return vals - + end """ $(SIGNATURES) @@ -462,54 +469,54 @@ performs the actual time-expansion in a sparse way. Returns SparseMatrixCSC """ -function time_expand(Xorg,term,tbl) +function time_expand(Xorg, term, tbl) # this is the predefined eventfield, usually "latency" tbl = DataFrame(tbl) onsets = Float64.(tbl[:, term.eventfields[1]])::Vector{Float64} # XXX if we have integer onsets, we could directly speed up matrix generation maybe? if typeof(term.eventfields) <: Array && length(term.eventfields) == 1 - bases = kernel.(Ref(term.basisfunction),onsets) + bases = kernel.(Ref(term.basisfunction), onsets) else - bases = kernel.(Ref(term.basisfunction),eachrow(tbl[!, term.eventfields])) + bases = kernel.(Ref(term.basisfunction), eachrow(tbl[!, term.eventfields])) end - return time_expand(Xorg,term,onsets, bases) + return time_expand(Xorg, term, onsets, bases) end function time_expand(Xorg, term, onsets, bases) - ncolsBasis = size(kernel(term.basisfunction,0), 2)::Int64 + ncolsBasis = size(kernel(term.basisfunction, 0), 2)::Int64 X = reshape(Xorg, size(Xorg, 1), :) # why is this necessary? ncolsX = size(X)[2]::Int64 - - - - - rows = timeexpand_rows(onsets,bases,shiftOnset(term.basisfunction),ncolsX) - cols = timeexpand_cols(term,bases,ncolsBasis,ncolsX) - vals = timeexpand_vals(bases,X,size(cols),ncolsX) - + + + + rows = timeexpand_rows(onsets, bases, shiftOnset(term.basisfunction), ncolsX) + cols = timeexpand_cols(term, bases, ncolsBasis, ncolsX) + + vals = timeexpand_vals(bases, X, size(cols), ncolsX) + #vals = vcat(vals...) ix = rows .> 0 #.&& vals .!= 0. A = @views sparse(rows[ix], cols[ix], vals[ix]) dropzeros!(A) - + return A end """ Helper function to decide whether all bases have the same number of columns per event """ -time_expand_allBasesSameCols(b::FIRBasis,bases,ncolBasis) = true # FIRBasis is always fast! -function time_expand_allBasesSameCols(basisfunction,bases,ncolsBasis) -fastpath = true -for b = eachindex(bases) - if length(unique(length.(nzrange.(Ref(bases[b]), 1:ncolsBasis)))) != 1 - return false +time_expand_allBasesSameCols(b::FIRBasis, bases, ncolBasis) = true # FIRBasis is always fast! +function time_expand_allBasesSameCols(basisfunction, bases, ncolsBasis) + fastpath = true + for b in eachindex(bases) + if length(unique(length.(nzrange.(Ref(bases[b]), 1:ncolsBasis)))) != 1 + return false + end end -end -return true + return true end """ @@ -558,14 +565,14 @@ end -function Base.show(io::IO,d::DesignMatrix) - println(io,"Unfold.DesignMatrix") - println(io,"Formulas: $(d.formulas)") +function Base.show(io::IO, d::DesignMatrix) + println(io, "Unfold.DesignMatrix") + println(io, "Formulas: $(d.formulas)") #display(io,d.formulas) - sz_evts = isa(d.events,Vector) ? size.(d.events) : size(d.events) - sz_Xs = (isa(d.Xs,Vector) | isa(d.Xs,Tuple)) ? size.(d.Xs) : size(d.Xs) - - println(io,"\nSizes: Xs: $sz_Xs, events: $sz_evts") - println(io,"\nuseful functions: formula(d),modelmatrix(d)") - println(io,"Fields: .formulas, .Xs, .events") -end \ No newline at end of file + sz_evts = isa(d.events, Vector) ? size.(d.events) : size(d.events) + sz_Xs = (isa(d.Xs, Vector) | isa(d.Xs, Tuple)) ? size.(d.Xs) : size(d.Xs) + + println(io, "\nSizes: Xs: $sz_Xs, events: $sz_evts") + println(io, "\nuseful functions: formula(d),modelmatrix(d)") + println(io, "Fields: .formulas, .Xs, .events") +end diff --git a/src/effects.jl b/src/effects.jl index 7bcc09e5..5c447349 100644 --- a/src/effects.jl +++ b/src/effects.jl @@ -1,8 +1,8 @@ import Effects: effects -import Effects:expand_grid -import Effects:typify +import Effects: expand_grid +import Effects: typify import Effects.typify -import Effects:_symequal +import Effects: _symequal import StatsModels.collect_matrix_terms import Base.getproperty using Effects @@ -23,98 +23,121 @@ Calculates marginal effects for all term-combinations in `design`. julia> effects(d,uf) ``` will result in 6 predicted values: A/-2, A/0, A/2, B/-2, B/0, B/2. -""" +""" -function effects(design::AbstractDict, model::UnfoldModel;typical=mean) +function effects(design::AbstractDict, model::UnfoldModel; typical = mean) reference_grid = expand_grid(design) form = Unfold.formula(model) # get formula # replace non-specified fields with "constants" - m = modelmatrix(model,false) # get the modelmatrix without timeexpansion - @debug "type form[1]",typeof(form[1]) - form_typical = _typify(reference_grid,form,m,typical) + m = modelmatrix(model, false) # get the modelmatrix without timeexpansion + @debug "type form[1]", typeof(form[1]) + form_typical = _typify(reference_grid, form, m, typical) @debug "type form_typical[1]", typeof(form_typical[1]) - eff = yhat(model,form_typical,reference_grid) - + eff = yhat(model, form_typical, reference_grid) + # because coefficients are 2D/3D arry, we have to cast it correctly to one big dataframe - if isa(eff,Tuple) + if isa(eff, Tuple) # TimeContinuous Model, we also get back other things like times & fromWhereToWhere a BasisFunction goes if typeof(form) <: FormulaTerm form = [form] end # figure out the basisname - bnames = [[form[ix].rhs.basisfunction.name] for ix in 1:length(form)] + bnames = [[form[ix].rhs.basisfunction.name] for ix = 1:length(form)] # repeat it as much as necessary - bnames = repeat.(bnames,[e.stop+e.step-1 for e in eff[1]]) + bnames = repeat.(bnames, [e.stop + e.step - 1 for e in eff[1]]) - result = DataFrame(cast_referenceGrid(reference_grid,eff[3],eff[2] ;basisname=vcat(bnames...))) - select!(result,Not(:latency)) # remove the latency column if it was added - elseif isa(eff,Vector) + result = DataFrame( + cast_referenceGrid(reference_grid, eff[3], eff[2]; basisname = vcat(bnames...)), + ) + select!(result, Not(:latency)) # remove the latency column if it was added + elseif isa(eff, Vector) # mass univariate with multiple effects - results_tmp = DataFrame.(cast_referenceGrid.(Ref(reference_grid),eff,Ref(times(model)))) + results_tmp = + DataFrame.(cast_referenceGrid.(Ref(reference_grid), eff, Ref(times(model)))) names = collect(keys(Unfold.design(model))) - [df.basisname .= n for (df,n) in zip(results_tmp,names)] + [df.basisname .= n for (df, n) in zip(results_tmp, names)] result = vcat(results_tmp...) else # normal mass univariate model - result = DataFrame(cast_referenceGrid(reference_grid,eff,times(model) )) + result = DataFrame(cast_referenceGrid(reference_grid, eff, times(model))) end - -return result + + return result end - - Effects.typify(reference_grid,form::AbstractArray,X;kwargs...) = typify.(Ref(reference_grid),form,Ref(X);kwargs...) - - # cast single form to a vector -_typify(reference_grid,form::FormulaTerm{<:InterceptTerm,<:Unfold.TimeExpandedTerm},m,typical) = _typify(reference_grid,[form],[m],typical) +Effects.typify(reference_grid, form::AbstractArray, X; kwargs...) = + typify.(Ref(reference_grid), form, Ref(X); kwargs...) + + +# cast single form to a vector +_typify( + reference_grid, + form::FormulaTerm{<:InterceptTerm,<:Unfold.TimeExpandedTerm}, + m, + typical, +) = _typify(reference_grid, [form], [m], typical) -function _typify(reference_grid, form::AbstractArray{<:FormulaTerm{<:InterceptTerm,<:Unfold.TimeExpandedTerm}},m::Vector{<:AbstractMatrix},typical) +function _typify( + reference_grid, + form::AbstractArray{<:FormulaTerm{<:InterceptTerm,<:Unfold.TimeExpandedTerm}}, + m::Vector{<:AbstractMatrix}, + typical, +) @debug "_typify - stripping away timeexpandedterm" - form_typical = Array{Any}(undef,1, length(form)) + form_typical = Array{Any}(undef, 1, length(form)) for f = 1:length(form) - + # strip of basisfunction and put it on afterwards again tmpf = deepcopy(form[f]) - + # create a Formula without Timeexpansion - tmpf = FormulaTerm(tmpf.lhs,tmpf.rhs.term) + tmpf = FormulaTerm(tmpf.lhs, tmpf.rhs.term) # typify that - tmpf = typify(reference_grid, tmpf, m[f]; typical=typical) - + tmpf = typify(reference_grid, tmpf, m[f]; typical = typical) + # regenerate TimeExpansion - tmpf = Unfold.TimeExpandedTerm(tmpf,form[f].rhs.basisfunction;eventfields=form[f].rhs.eventfields) + tmpf = Unfold.TimeExpandedTerm( + tmpf, + form[f].rhs.basisfunction; + eventfields = form[f].rhs.eventfields, + ) form_typical[f] = tmpf end return form_typical - end - function _typify(reference_grid,form::FormulaTerm,m,typical) - - return [typify(reference_grid, form, m; typical=typical)] +end +function _typify(reference_grid, form::FormulaTerm, m, typical) + + return [typify(reference_grid, form, m; typical = typical)] - end +end - function _typify(reference_grid,form::AbstractArray{<:FormulaTerm},m::Vector{<:AbstractMatrix},typical) +function _typify( + reference_grid, + form::AbstractArray{<:FormulaTerm}, + m::Vector{<:AbstractMatrix}, + typical, +) # Mass Univariate with multiple effects @debug "_typify going the mass univariate route" out = [] for k = 1:length(form) - push!(out,typify(reference_grid, form[k], m[k]; typical=typical)) + push!(out, typify(reference_grid, form[k], m[k]; typical = typical)) end return out - end - -function cast_referenceGrid(r,eff,times;basisname=nothing) +end + +function cast_referenceGrid(r, eff, times; basisname = nothing) nchan = size(eff, 2) # correct - neff = size(r,1) # how many effects requested - neffCol = size(r,2) # how many predictors - ncols = size(eff,1) ÷ neff # typically ntimes + neff = size(r, 1) # how many effects requested + neffCol = size(r, 2) # how many predictors + ncols = size(eff, 1) ÷ neff # typically ntimes + + - - # replicate # for each predictor in r (reference grid), we need this at the bottom @@ -123,30 +146,35 @@ function cast_referenceGrid(r,eff,times;basisname=nothing) else nbases = length(unique(basisname)) end - coefs_rep = Array{Array}(undef,nbases,neffCol) - - + coefs_rep = Array{Array}(undef, nbases, neffCol) + + for k = 1:neffCol # in case we have only a single basis (e.g. mass univariate), we can directly fill in all values ixList = [] if isnothing(basisname) - ix = ones(ncols) .== 1. - append!(ixList,[ix]) + ix = ones(ncols) .== 1.0 + append!(ixList, [ix]) else #in case of multiple bases, we have to do it iteratively, because the bases can be different length - for b = unique(basisname) - ix = basisname[1:neff:end].==b - append!(ixList,[ix]) + for b in unique(basisname) + ix = basisname[1:neff:end] .== b + append!(ixList, [ix]) end end for i_ix = 1:length(ixList) - - coefs_rep[i_ix,k] = linearize(permutedims(repeat(r[:,k], outer = [1, nchan, sum(ixList[i_ix])]), [2,3,1])) + + coefs_rep[i_ix, k] = linearize( + permutedims( + repeat(r[:, k], outer = [1, nchan, sum(ixList[i_ix])]), + [2, 3, 1], + ), + ) end end - + # often the "times" vector - if length(times) == neff*ncols + if length(times) == neff * ncols # in case we have timeexpanded, times is already in long format and doesnt need to be repeated for each coefficient colnames_basis_rep = permutedims(repeat(times, 1, nchan, 1), [2 1 3]) else @@ -155,23 +183,25 @@ function cast_referenceGrid(r,eff,times;basisname=nothing) # for multiple channels chan_rep = repeat(1:nchan, 1, ncols, neff) - + # for mass univariate there is no basisname if isnothing(basisname) - basisname = fill(nothing,ncols) - basisname_rep = permutedims(repeat(basisname,1,nchan,neff),[2,1,3]) + basisname = fill(nothing, ncols) + basisname_rep = permutedims(repeat(basisname, 1, nchan, neff), [2, 1, 3]) else - basisname_rep = permutedims(repeat(basisname,1,nchan,1),[2,1,3]) + basisname_rep = permutedims(repeat(basisname, 1, nchan, 1), [2, 1, 3]) end - - result = Dict( :yhat => linearize(eff'), - :time => linearize(colnames_basis_rep), - :channel => linearize(chan_rep), - :basisname =>linearize(basisname_rep)) + + result = Dict( + :yhat => linearize(eff'), + :time => linearize(colnames_basis_rep), + :channel => linearize(chan_rep), + :basisname => linearize(basisname_rep), + ) for k = 1:neffCol - push!(result,Symbol(names(r)[k]) =>vcat(vcat(coefs_rep[:,k]...)...)) + push!(result, Symbol(names(r)[k]) => vcat(vcat(coefs_rep[:, k]...)...)) end return result @@ -179,7 +209,7 @@ end #Effects._symequal(t1::AbstractTerm,t2::Unfold.TimeExpandedTerm) = _symequal(t1,t2.term) #function Effects._replace(matrix_term::MatrixTerm{<:Tuple{<:Unfold.TimeExpandedTerm}},typicals::Dict) - + # replaced_term = MatrixTerm((Effects._replace.(matrix_term.terms, Ref(typicals))...,)) # basisfunctionTerm = getfield(matrix_term,:terms)[1] # return TimeExpandedTerm(replaced_term,basisfunctionTerm.basisfunction,basisfunctionTerm.eventfields) diff --git a/src/eventhandling.jl b/src/eventhandling.jl index dcc13027..71bf7642 100644 --- a/src/eventhandling.jl +++ b/src/eventhandling.jl @@ -8,46 +8,47 @@ copy field-info from source to closest target julia> # copy reaction time values from button press to closest stimulus immediately before button press julia> copy_eventinfo!(evts,"button"=>"stimulus",:reaction_time;match_fun="s after t") """ -function copy_eventinfo!(evts,fromTo,field;match_fun="closest") +function copy_eventinfo!(evts, fromTo, field; match_fun = "closest") source = fromTo.first target = fromTo.second - if isa(field,Pair) + if isa(field, Pair) source_field = field.first target_field = field.second else source_field = field target_field = field end - source_ix = findall(isequal(source),evts.trial_type) - target_ix = findall(isequal(target),evts.trial_type) + source_ix = findall(isequal(source), evts.trial_type) + target_ix = findall(isequal(target), evts.trial_type) - isempty(source_ix) && error("couldnt find source entries ($source) in evts.trial_type" ) - isempty(target_ix) && error("couldnt find target entries ($target) in evts.trial_type" ) + isempty(source_ix) && error("couldnt find source entries ($source) in evts.trial_type") + isempty(target_ix) && error("couldnt find target entries ($target) in evts.trial_type") # for some matching functions we want to find the minimum, but no negative number - filter_greaterzero = x-> x >= 0 ? x : Inf - + filter_greaterzero = x -> x >= 0 ? x : Inf + if match_fun == "closest" - match_fun = (a,b)-> argmin(abs.(a .- b)) # find closest before or after + match_fun = (a, b) -> argmin(abs.(a .- b)) # find closest before or after elseif match_fun == "s after t" || match_fun == "t before s" - match_fun = (s,t)-> argmin(filter_greaterzero.(s .- t)) # source after target + match_fun = (s, t) -> argmin(filter_greaterzero.(s .- t)) # source after target elseif match_fun == "s before t" || match_fun == "t after s" - match_fun = (s,t)-> argmin(filter_greaterzero.(t .- s)) # source before target + match_fun = (s, t) -> argmin(filter_greaterzero.(t .- s)) # source before target elseif !callable(match_fun) error("bad match_fun input") end - ix = [match_fun(evts[s,:latency], evts[target_ix,:latency]) for s in source_ix] + ix = [match_fun(evts[s, :latency], evts[target_ix, :latency]) for s in source_ix] ix = target_ix[ix] - payload = evts[source_ix,source_field] + payload = evts[source_ix, source_field] # if targetfield does not exist, create one with Union(missing,target_field_datatype) if !any(names(evts) .== target_field) - evts[!,target_field] = Array{Union{Missing,typeof(payload[1])}}(missing,nrow(evts)) + evts[!, target_field] = + Array{Union{Missing,typeof(payload[1])}}(missing, nrow(evts)) end - + # fill it in - evts[ix,target_field] .= payload + evts[ix, target_field] .= payload return evts -end \ No newline at end of file +end diff --git a/src/fit.jl b/src/fit.jl index 298c7d4b..f7d420b2 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -43,10 +43,10 @@ function StatsModels.fit( tbl::DataFrame, data::AbstractArray; kwargs..., - ) +) detectedType = designToModeltype(design) - fit(detectedType,design,tbl,data;kwargs...) + fit(detectedType, design, tbl, data; kwargs...) end @@ -57,7 +57,7 @@ function StatsModels.fit( data::AbstractArray; kwargs..., ) -fit(UnfoldModelType(design),design,tbl,data;kwargs...) + fit(UnfoldModelType(design), design, tbl, data; kwargs...) end @@ -67,8 +67,8 @@ function StatsModels.fit( tbl::DataFrame, data::AbstractArray; kwargs..., - ) - +) + designmatrix!(uf, tbl; kwargs...) fit!(uf, data; kwargs...) @@ -100,8 +100,8 @@ isMixedModelFormula(f::InteractionTerm) = false isMixedModelFormula(f::ConstantTerm) = false isMixedModelFormula(f::Term) = false #isMixedModelFormula(f::FunctionTerm) = false -function isMixedModelFormula(f::FunctionTerm) - try +function isMixedModelFormula(f::FunctionTerm) + try isMixedModelFormula(f.f) catch isMixedModelFormula(f.forig) # StatsMoels <0.7 @@ -157,7 +157,8 @@ function StatsModels.fit( tbl::DataFrame, data::AbstractMatrix, args...; - kwargs...)where {T<:Union{UnfoldLinearMixedModel,UnfoldLinearModel}} + kwargs..., +) where {T<:Union{UnfoldLinearMixedModel,UnfoldLinearModel}} @debug("MassUnivariate data array is size (X,Y), reshaping to (1,X,Y)") data = reshape(data, 1, size(data)...) return fit(ufmodel, design, tbl, data, args...; kwargs...) @@ -175,49 +176,48 @@ function StatsModels.fit!( @assert ~isempty(designmatrix(uf)) @assert typeof(first(values(design(uf)))[1]) <: FormulaTerm "InputError in design(uf) - :key=>(FORMULA,basis/times), formula not found. Maybe formula wasn't at the first place?" - @assert (typeof(first(values(design(uf)))[2]) <: AbstractVector) ⊻ (typeof(uf) <: UnfoldLinearModelContinuousTime) "InputError: Either a basis function was declared, but a UnfoldLinearModel was built, or a times-vector (and no basis function) was given, but a UnfoldLinearModelContinuousTime was asked for." - if isa(uf,UnfoldLinearModel) - @assert length(first(values(design(uf)))[2]) == size(data,length(size(data))-1) "Times Vector does not match second last dimension of input data - forgot to epoch?" + @assert (typeof(first(values(design(uf)))[2]) <: AbstractVector) ⊻ + (typeof(uf) <: UnfoldLinearModelContinuousTime) "InputError: Either a basis function was declared, but a UnfoldLinearModel was built, or a times-vector (and no basis function) was given, but a UnfoldLinearModelContinuousTime was asked for." + if isa(uf, UnfoldLinearModel) + @assert length(first(values(design(uf)))[2]) == size(data, length(size(data)) - 1) "Times Vector does not match second last dimension of input data - forgot to epoch?" end - + X = modelmatrix(uf) @debug "UnfoldLinearModel(ContinuousTime), datasize: $(size(data))" - - if isa(uf,UnfoldLinearModel) + + if isa(uf, UnfoldLinearModel) d = designmatrix(uf) - if isa(X,Vector) - # mass univariate with multiple events fitted at the same time - - coefs = [] - for m = 1:length(X) - # the main issue is, that the designmatrices are subsets of the event table - we have - # to do the same for the data, but data and designmatrix dont know much about each other. - # Thus we use parentindices() to get the original indices of the @view events[...] from desigmatrix.jl - push!(coefs,solver(X[m], @view data[:,:,parentindices(d.events[m])[1]])) - end - @debug @show [size(c.estimate) for c in coefs] - uf.modelfit = LinearModelFit( - cat([c.estimate for c in coefs]...,dims=3), - [c.info for c in coefs], - cat([c.standarderror for c in coefs]...,dims=3) - ) - return # we are done here - - elseif isa(d.events,SubDataFrame) + if isa(X, Vector) + # mass univariate with multiple events fitted at the same time + + coefs = [] + for m = 1:length(X) + # the main issue is, that the designmatrices are subsets of the event table - we have + # to do the same for the data, but data and designmatrix dont know much about each other. + # Thus we use parentindices() to get the original indices of the @view events[...] from desigmatrix.jl + push!(coefs, solver(X[m], @view data[:, :, parentindices(d.events[m])[1]])) + end + @debug @show [size(c.estimate) for c in coefs] + uf.modelfit = LinearModelFit( + cat([c.estimate for c in coefs]..., dims = 3), + [c.info for c in coefs], + cat([c.standarderror for c in coefs]..., dims = 3), + ) + return # we are done here + + elseif isa(d.events, SubDataFrame) # in case the user specified an event to subset (and not any) we have to use the view from now on - data = @view data[:,:,parentindices(d.events)[1]] + data = @view data[:, :, parentindices(d.events)[1]] end end - # mass univariate, data = ch x times x epochs - X, data = zeropad(X, data) + # mass univariate, data = ch x times x epochs + X, data = zeropad(X, data) - uf.modelfit = solver(X, data) - return + uf.modelfit = solver(X, data) + return end - - diff --git a/src/io.jl b/src/io.jl index 6df04b03..16bb51df 100644 --- a/src/io.jl +++ b/src/io.jl @@ -26,9 +26,17 @@ Save UnfoldModel in a (by default uncompressed) .jld2 file. For memory efficiency the designmatrix is set to missing. If needed, it can be reconstructed when loading the model. """ -function FileIO.save(file, uf::T; compress=false) where {T<:UnfoldModel} - jldopen(file, "w"; compress=compress) do f - f["uf"] = T(uf.design, Unfold.DesignMatrix(designmatrix(uf).formulas, missing, designmatrix(uf).events), uf.modelfit) +function FileIO.save(file, uf::T; compress = false) where {T<:UnfoldModel} + jldopen(file, "w"; compress = compress) do f + f["uf"] = T( + uf.design, + Unfold.DesignMatrix( + designmatrix(uf).formulas, + missing, + designmatrix(uf).events, + ), + uf.modelfit, + ) end end @@ -41,7 +49,7 @@ Load UnfoldModel from a .jld2 file. By default, the designmatrix is reconstructed. If it is not needed set `generate_Xs=false` which improves time-efficiency. """ -function FileIO.load(file, ::Type{<:UnfoldModel}; generate_Xs=true) +function FileIO.load(file, ::Type{<:UnfoldModel}; generate_Xs = true) f = jldopen(file, "r") uf = f["uf"] close(f) @@ -58,4 +66,4 @@ function FileIO.load(file, ::Type{<:UnfoldModel}; generate_Xs=true) # reintegrate the designmatrix return typeof(uf)(uf.design, Unfold.DesignMatrix(form, X, events), uf.modelfit) -end \ No newline at end of file +end diff --git a/src/predict.jl b/src/predict.jl index cd851e79..38fc3a76 100644 --- a/src/predict.jl +++ b/src/predict.jl @@ -6,22 +6,22 @@ using DocStringExtensions: format function StatsBase.predict(model::UnfoldModel, events) # make a copy of it so we don't change it outside the function newevents = copy(events) - + formulas = formula(model) if typeof(formulas) <: FormulaTerm formulas = [formulas] end if typeof(model) == UnfoldLinearModel - eff = yhat(model,formulas[1],newevents) - timesVec = gen_timeev(times(model),size(newevents, 1)) + eff = yhat(model, formulas[1], newevents) + timesVec = gen_timeev(times(model), size(newevents, 1)) else - fromTo,timesVec,eff = yhat(model,formulas,newevents) + fromTo, timesVec, eff = yhat(model, formulas, newevents) end - + # init the meta dataframe metaData = DataFrame([:time => vcat(timesVec...), :basisname => ""]) - + for c in names(newevents) metaData[:, c] .= newevents[1, c] # assign first element in order to have same column type end @@ -41,31 +41,31 @@ function StatsBase.predict(model::UnfoldModel, events) # shift variable to keep track of multiple basisfunctions shift = 0 # for each basis function - + for (bIx, basisfun) in enumerate(fromTo) - + # go through all predictors for (i, fstart) in enumerate(basisfun[1:end]) - - fend = basisfun[i] + basisfun.step-1 - + + fend = basisfun[i] + basisfun.step - 1 + # couldn't figure out how to broadcast everything directly (i.e out[fstart:fend,names(newevents)] .= newevents[i,:]) # copy the correct metadata - + for j = fstart:fend metaData[j+shift, names(newevents)] = newevents[i, :] end # add basisfunction name metaData[shift.+(fstart:fend), :basisname] .= formulas[bIx].rhs.basisfunction.name - + end # the next meta data has to be at the end - shift += basisfun[end]-1+basisfun.step + shift += basisfun[end] - 1 + basisfun.step end end - + out = DataFrame([:yhat => vec(reshape(eff, :, 1))]) nchannel = size(eff, 2) @@ -87,16 +87,16 @@ function yhat( X = modelcols.(formulas, Ref(newevents)) co = coef(model) - Xsizes = size.(X,Ref(2)) - Xsizes_cumsum = vcat(0,cumsum(Xsizes)) - + Xsizes = size.(X, Ref(2)) + Xsizes_cumsum = vcat(0, cumsum(Xsizes)) + indexes = [(Xsizes_cumsum[ix]+1):Xsizes_cumsum[ix+1] for ix = 1:length(formulas)] - - coArray = [co[:,:,ix] for ix in indexes] - - return yhat.(coArray,X) - + coArray = [co[:, :, ix] for ix in indexes] + + + return yhat.(coArray, X) + end @@ -116,24 +116,25 @@ function yhat( end -yhat(model::UnfoldLinearModelContinuousTime,formulas::MatrixTerm,events) = yhat(model,formulas.terms,events) +yhat(model::UnfoldLinearModelContinuousTime, formulas::MatrixTerm, events) = + yhat(model, formulas.terms, events) -function yhat(model::UnfoldLinearModelContinuousTime,formulas,events)#::AbstractArray{AbstractTerm}) +function yhat(model::UnfoldLinearModelContinuousTime, formulas, events)#::AbstractArray{AbstractTerm}) X = [] fromTo = [] timesVec = [] for f in formulas # due to many reasons, there can be several different options of how formula arrives here - not easy to catch via multiple dispatch - if !(isa(f,TimeExpandedTerm)) - if !(isa(f,MatrixTerm)) - f = f.rhs - #elseif !(isa(f,Effects.TypicalTerm)) - + if !(isa(f, TimeExpandedTerm)) + if !(isa(f, MatrixTerm)) + f = f.rhs + #elseif !(isa(f,Effects.TypicalTerm)) + else f = f.terms[1] end - + end # find out how long each designmatrix is n_range = length(times(f.basisfunction)) @@ -144,10 +145,10 @@ function yhat(model::UnfoldLinearModelContinuousTime,formulas,events)#::Abstract range(-n_negative + 1, step = n_range, length = size(events, 1)) - + # get the model Xsingle = modelcols(f, events) - + timesSingle = times(f.basisfunction) # this pertains only to FIR-models @@ -161,18 +162,22 @@ function yhat(model::UnfoldLinearModelContinuousTime,formulas,events)#::Abstract # f&g are between two samples, thus the design matrix would interpolate between them. Thus has as a result, that the designmatrix is +1 longer than what would naively be expected # # because in "predict" we define where samples onset, we can remove the last sample, it s always 0 anyway, but to be sure we test it - + if typeof(f.basisfunction) <: FIRBasis - - keep = ones(size(Xsingle,1)) - keep[range(length(timesSingle), size(Xsingle,1),step=length(timesSingle))] .= 0 + + keep = ones(size(Xsingle, 1)) + keep[range( + length(timesSingle), + size(Xsingle, 1), + step = length(timesSingle), + )] .= 0 Xsingle = Xsingle[keep.==1, :] timesSingle = timesSingle[1:end-1] - n_range = n_range-1 + n_range = n_range - 1 end # combine designmats append!(X, [Xsingle]) - + # keep track of how long each event is append!(fromTo, [range(1, step = n_range, length = size(events, 1))]) # keep track of the times @@ -184,14 +189,18 @@ function yhat(model::UnfoldLinearModelContinuousTime,formulas,events)#::Abstract Xconcat = blockdiag(X...) # calculate yhat - eff = yhat(model,Xconcat) + eff = yhat(model, Xconcat) # output is a bit ugly, but we need the other two vectors as well. Should be refactored at some point - but not right now ;-) #XXX - return fromTo,timesVec,eff - + return fromTo, timesVec, eff + end -function yhat(model::UnfoldLinearModelContinuousTime,X::AbstractArray{T,2};times=nothing) where {T<:Union{Missing, <:Number}} - +function yhat( + model::UnfoldLinearModelContinuousTime, + X::AbstractArray{T,2}; + times = nothing, +) where {T<:Union{Missing,<:Number}} + yhat = X * coef(model)' return yhat @@ -213,16 +222,11 @@ function yhat( # function that calculates coef*designmat, but in the ch x times x coef vector # setup the output matrix, has to be a matrix # then transforms it back to 2D matrix times/coef x ch to be compatible with the timecontinuous format - yhat = Array{Union{Missing,Float64}}( - missing, - size(coef, 1), - size(X, 1), - size(coef, 2), - ) + yhat = Array{Union{Missing,Float64}}(missing, size(coef, 1), size(X, 1), size(coef, 2)) for ch = 1:size(coef, 1) yhat[ch, :, :] = X * permutedims(coef[ch, :, :], (2, 1)) end - + # bring the yhat into a ch x yhat format yhat = reshape(permutedims(yhat, (1, 3, 2)), size(yhat, 1), :) yhat = permutedims(yhat, (2, 1)) @@ -236,7 +240,7 @@ times( ) = first(values(d))[2]#[k[2] for k in values(d)] # probably going for steprange would be better -function gen_timeev(timesVec,nRows) +function gen_timeev(timesVec, nRows) timesVec = repeat(timesVec, nRows) return timesVec -end \ No newline at end of file +end diff --git a/src/solver.jl b/src/solver.jl index cc32c048..e7faae97 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -2,22 +2,23 @@ using StatsBase: var function solver_default( X, data::AbstractArray{T,2}; - stderror=false, - multithreading=true, - showprogress=true, + stderror = false, + multithreading = true, + showprogress = true, ) where {T<:Union{Missing,<:Number}} minfo = Array{IterativeSolvers.ConvergenceHistory,1}(undef, size(data, 1)) beta = zeros(size(data, 1), size(X, 2)) # had issues with undef - p = Progress(size(data, 1); enabled=showprogress) + p = Progress(size(data, 1); enabled = showprogress) @maybe_threads multithreading for ch = 1:size(data, 1) dd = view(data, ch, :) ix = @. !ismissing(dd) # use the previous channel as a starting point ch == 1 || copyto!(view(beta, ch, :), view(beta, ch - 1, :)) - beta[ch, :], h = lsmr!(@view(beta[ch, :]), (X[ix, :]), @view(data[ch, ix]), log=true) + beta[ch, :], h = + lsmr!(@view(beta[ch, :]), (X[ix, :]), @view(data[ch, ix]), log = true) minfo[ch] = h next!(p) @@ -36,7 +37,7 @@ end function calculate_stderror(Xdc, data::Matrix{T}, beta) where {T<:Union{Missing,<:Number}} # remove missings - ix = any(.!ismissing.(data), dims=1)[1, :] + ix = any(.!ismissing.(data), dims = 1)[1, :] if length(ix) != size(data, 2) @warn( "Limitation: Missing data are calculated over all channels for standard error" @@ -93,13 +94,13 @@ end function solver_default( X, data::AbstractArray{T,3}; - stderror=false, - multithreading=true, - showprogress=true, + stderror = false, + multithreading = true, + showprogress = true, ) where {T<:Union{Missing,<:Number}} #beta = Array{Union{Missing,Number}}(undef, size(data, 1), size(data, 2), size(X, 2)) beta = zeros(Union{Missing,Number}, size(data, 1), size(data, 2), size(X, 2)) - p = Progress(size(data, 1); enabled=showprogress) + p = Progress(size(data, 1); enabled = showprogress) @maybe_threads multithreading for ch = 1:size(data, 1) for t = 1:size(data, 2) # @debug("$(ndims(data,)),$t,$ch") @@ -122,13 +123,13 @@ function solver_default( return modelfit end -solver_b2b(X, data, cross_val_reps) = solver_b2b(X, data, cross_val_reps=cross_val_reps) +solver_b2b(X, data, cross_val_reps) = solver_b2b(X, data, cross_val_reps = cross_val_reps) function solver_b2b( X, data::AbstractArray{T,3}; - cross_val_reps=10, - multithreading=true, - showprogress=true, + cross_val_reps = 10, + multithreading = true, + showprogress = true, ) where {T<:Union{Missing,<:Number}} X, data = dropMissingEpochs(X, data) @@ -137,7 +138,7 @@ function solver_b2b( E = zeros(size(data, 2), size(X, 2), size(X, 2)) W = Array{Float64}(undef, size(data, 2), size(X, 2), size(data, 1)) - prog = Progress(size(data, 2) * cross_val_reps, 0.1; enabled=showprogress) + prog = Progress(size(data, 2) * cross_val_reps, 0.1; enabled = showprogress) @maybe_threads multithreading for m = 1:cross_val_reps k_ix = collect(Kfold(size(data, 3), 2)) X1 = @view X[k_ix[1], :] @@ -153,7 +154,7 @@ function solver_b2b( H = X2 \ (Y2' * G) E[t, :, :] += Diagonal(H[diagind(H)]) - ProgressMeter.next!(prog; showvalues=[(:time, t), (:cross_val_rep, m)]) + ProgressMeter.next!(prog; showvalues = [(:time, t), (:cross_val_rep, m)]) end E[t, :, :] = E[t, :, :] ./ cross_val_reps W[t, :, :] = (X * E[t, :, :])' / data[:, t, :] @@ -161,10 +162,9 @@ function solver_b2b( end # extract diagonal - beta = mapslices(diag, E, dims=[2, 3]) + beta = mapslices(diag, E, dims = [2, 3]) # reshape to conform to ch x time x pred beta = permutedims(beta, [3 1 2]) modelinfo = Dict("W" => W, "E" => E, "cross_val_reps" => cross_val_reps) # no history implemented (yet?) return LinearModelFit(beta, modelinfo) end - diff --git a/src/splinepredictors.jl b/src/splinepredictors.jl index e69de29b..8b137891 100644 --- a/src/splinepredictors.jl +++ b/src/splinepredictors.jl @@ -0,0 +1 @@ + diff --git a/src/statistics.jl b/src/statistics.jl index e69de29b..8b137891 100644 --- a/src/statistics.jl +++ b/src/statistics.jl @@ -0,0 +1 @@ + diff --git a/src/timeexpandedterm.jl b/src/timeexpandedterm.jl index 2c58c5e8..b5949aa1 100644 --- a/src/timeexpandedterm.jl +++ b/src/timeexpandedterm.jl @@ -9,7 +9,7 @@ $(FIELDS) julia> b = TimeExpandedTerm(term,kernel,[:latencyTR,:durationTR]) ``` """ -struct TimeExpandedTerm{T<:AbstractTerm} <: AbstractTerm +struct TimeExpandedTerm{T<:AbstractTerm} <: AbstractTerm "Term that the basis function is applied to. This is regularly called in other functions to get e.g. term-coefnames and timeexpand those" term::T "Kernel that determines what should happen to the designmatrix of the term" @@ -39,8 +39,10 @@ StatsModels.width(term::TimeExpandedTerm) = width(term.basisfunction) StatsModels.terms(t::TimeExpandedTerm) = terms(t.term) function Base.show(io::IO, p::TimeExpandedTerm) - print(io, "$(p.basisfunction.name): timeexpand($(p.term)) for times $(times(p.basisfunction))") - + print( + io, + "$(p.basisfunction.name): timeexpand($(p.term)) for times $(times(p.basisfunction))", + ) + #println(io, "$(coefnames(p))") end - diff --git a/src/typedefinitions.jl b/src/typedefinitions.jl index bbedd0c8..3b6afbd5 100644 --- a/src/typedefinitions.jl +++ b/src/typedefinitions.jl @@ -82,4 +82,4 @@ function Base.show(io::IO, obj::UnfoldModel) end -abstract type AbstractSplineTerm <:AbstractTerm end +abstract type AbstractSplineTerm <: AbstractTerm end diff --git a/src/utilities.jl b/src/utilities.jl index cbae603e..01768cad 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -12,19 +12,13 @@ function epoch( kwargs..., ) """ -function epoch( - data::Array{T,1}, - tbl, - τ, - sfreq; - kwargs..., -) where {T<:Union{Missing,Number}} +function epoch(data::Array{T,1}, tbl, τ, sfreq; kwargs...) where {T<:Union{Missing,Number}} data_r = reshape(data, (1, :)) epoch(data_r, tbl, τ, sfreq; kwargs...) end function epoch(data::Matrix, tbl::DataFrame, τ::Vector, sfreq; kwargs...) - return epoch(data, tbl, (τ[1],τ[2]), sfreq; kwargs...) + return epoch(data, tbl, (τ[1], τ[2]), sfreq; kwargs...) end function epoch( @@ -33,8 +27,8 @@ function epoch( τ::Tuple{Number,Number}, sfreq; eventtime::String = "latency", -) -return epoch(data, tbl, τ, sfreq; eventtime = Symbol(eventtime)) +) + return epoch(data, tbl, τ, sfreq; eventtime = Symbol(eventtime)) end function epoch( @@ -68,9 +62,9 @@ function epoch( for si = 1:size(tbl, 1) #eventonset = tbl[si,eventtime] # in samples #d_start = eventonset - d_start = Int(round(tbl[si, eventtime]) + times[1].*sfreq) - d_end = Int(round(tbl[si, eventtime]) +times[end].*sfreq) - + d_start = Int(round(tbl[si, eventtime]) + times[1] .* sfreq) + d_end = Int(round(tbl[si, eventtime]) + times[end] .* sfreq) + e_start = 1 e_end = lenEpochs #println("d: $(size(data)),e: $(size(epochs)) | $d_start,$d_end,$e_start,$e_end | $(tbl[si,eventtime])") @@ -122,7 +116,10 @@ Equates the length of data and designmatrix by cutting the shorter one The reason we need this is because when generating the designmatrix, we do not know how long the data actually are. We only assume that event-latencies are synchronized with the data """ -function zeropad(X::AbstractMatrix, data::AbstractArray{T,2}) where {T<:Union{Missing,<:Number}} +function zeropad( + X::AbstractMatrix, + data::AbstractArray{T,2}, +) where {T<:Union{Missing,<:Number}} @debug("2d zeropad") if size(X, 1) > size(data, 2) X = X[1:size(data, 2), :] @@ -131,21 +128,27 @@ function zeropad(X::AbstractMatrix, data::AbstractArray{T,2}) where {T<:Union{Mi end return X, data end -function zeropad(X::AbstractMatrix, data::AbstractVector{T}) where {T<:Union{Missing,<:Number}} +function zeropad( + X::AbstractMatrix, + data::AbstractVector{T}, +) where {T<:Union{Missing,<:Number}} @debug("1d zeropad") if size(X, 1) > length(data) - X = X[1:length(data),:] + X = X[1:length(data), :] else data = data[1:size(X, 1)] end return X, data end -function zeropad(X::AbstractMatrix, data::AbstractArray{T,3}) where {T<:Union{Missing,<:Number}} +function zeropad( + X::AbstractMatrix, + data::AbstractArray{T,3}, +) where {T<:Union{Missing,<:Number}} @debug("3d zeropad") - + @assert size(X, 1) == size(data, 3) "Your events are not of the same size as your last dimension of data" - + return X, data end @@ -163,7 +166,7 @@ function clean_data( end -macro maybe_threads(multithreading,code) +macro maybe_threads(multithreading, code) return esc(:( if multithreading Threads.@threads($code) @@ -171,4 +174,4 @@ macro maybe_threads(multithreading,code) $code end )) -end \ No newline at end of file +end diff --git a/test/basisfunctions.jl b/test/basisfunctions.jl index 785b8307..44ccae9f 100644 --- a/test/basisfunctions.jl +++ b/test/basisfunctions.jl @@ -1,39 +1,39 @@ @testset "FIR" begin -firbase = firbasis((-1, 1), 10) + firbase = firbasis((-1, 1), 10) -# test optional call -@test Unfold.kernel(firbase,1) == Unfold.kernel(firbasis((-1, 1), 10),1) + # test optional call + @test Unfold.kernel(firbase, 1) == Unfold.kernel(firbasis((-1, 1), 10), 1) -@test typeof(Unfold.name(firbase)) <: String -# test basics of basisfunction -@test length(collect(Unfold.colnames(firbase))) == 21 -@test unique(Unfold.kernel(firbase,1)) == [1.0, 0.0] + @test typeof(Unfold.name(firbase)) <: String + # test basics of basisfunction + @test length(collect(Unfold.colnames(firbase))) == 21 + @test unique(Unfold.kernel(firbase, 1)) == [1.0, 0.0] -# test length consistency -@test length(Unfold.colnames(firbase)) ==size(Unfold.kernel(firbase,3.1))[2] -@test length(Unfold.times(firbase)) == size(Unfold.kernel(firbase,3.1))[1] + # test length consistency + @test length(Unfold.colnames(firbase)) == size(Unfold.kernel(firbase, 3.1))[2] + @test length(Unfold.times(firbase)) == size(Unfold.kernel(firbase, 3.1))[1] -# testing the non-sampling rate samples -@test Unfold.kernel(firbase,0.5)[1:3, 1:3] == [0.5 0.0 0.0; 0.5 0.5 0.0; 0.0 0.5 0.5] + # testing the non-sampling rate samples + @test Unfold.kernel(firbase, 0.5)[1:3, 1:3] == [0.5 0.0 0.0; 0.5 0.5 0.0; 0.0 0.5 0.5] end @testset "BOLD" begin -boldbase = hrfbasis(2.0, name = "test") + boldbase = hrfbasis(2.0, name = "test") -@test Unfold.name(boldbase) == "test" -@test Unfold.kernel(boldbase,0) == Unfold.kernel(boldbase,1) -@test Unfold.kernel(boldbase,0.1) != Unfold.kernel(boldbase,1) # testing fractional kernel generation -@test findmax(Unfold.kernel(boldbase,0.3))[2][1] == 4 + @test Unfold.name(boldbase) == "test" + @test Unfold.kernel(boldbase, 0) == Unfold.kernel(boldbase, 1) + @test Unfold.kernel(boldbase, 0.1) != Unfold.kernel(boldbase, 1) # testing fractional kernel generation + @test findmax(Unfold.kernel(boldbase, 0.3))[2][1] == 4 end @testset "timespline" begin - splinebase = Unfold.splinebasis(τ=(-1,1),sfreq=20,nsplines=10,name="basisA") - - @test length(Unfold.colnames(splinebase)) ==size(Unfold.kernel(splinebase,3.1))[2] - @test length(Unfold.times(splinebase)) == size(Unfold.kernel(splinebase,3.1))[1] + splinebase = Unfold.splinebasis(τ = (-1, 1), sfreq = 20, nsplines = 10, name = "basisA") -end \ No newline at end of file + @test length(Unfold.colnames(splinebase)) == size(Unfold.kernel(splinebase, 3.1))[2] + @test length(Unfold.times(splinebase)) == size(Unfold.kernel(splinebase, 3.1))[1] + +end diff --git a/test/designmatrix.jl b/test/designmatrix.jl index f6e2d60b..e2668072 100644 --- a/test/designmatrix.jl +++ b/test/designmatrix.jl @@ -16,227 +16,244 @@ shouldBePos[4, :] = [0, 0, 1, 0] @testset "basic designmat" begin -## test negative -basisfunction = firbasis(τ = (-3, 0), sfreq = 1, name = "testing") -timeexpandterm = Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, :latency); -Xdc = Unfold.time_expand(X, timeexpandterm, tbl) -@test all(isapprox.(Matrix(Xdc)[1:4, 1:4], shouldBeNeg, atol = 1e-15)) - -## Test Positive only -basisfunction = firbasis(τ = (1, 4), sfreq = 1, name = "testing") -timeexpandterm = Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, :latency); -Xdc = Unfold.time_expand(X, timeexpandterm, tbl) -println(Matrix(Xdc)) - -@test all(isapprox.(Matrix(Xdc)[1:4, 1:4], shouldBePos, atol = 1e-15)) + ## test negative + basisfunction = firbasis(τ = (-3, 0), sfreq = 1, name = "testing") + timeexpandterm = + Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, :latency) + Xdc = Unfold.time_expand(X, timeexpandterm, tbl) + @test all(isapprox.(Matrix(Xdc)[1:4, 1:4], shouldBeNeg, atol = 1e-15)) + + ## Test Positive only + basisfunction = firbasis(τ = (1, 4), sfreq = 1, name = "testing") + timeexpandterm = + Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, :latency) + Xdc = Unfold.time_expand(X, timeexpandterm, tbl) + println(Matrix(Xdc)) + + @test all(isapprox.(Matrix(Xdc)[1:4, 1:4], shouldBePos, atol = 1e-15)) end @testset "customized eventfields" begin -tbl2 = tbl = DataFrame([1 4]', [:onset]) - -timeexpandterm_latency = Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction); -timeexpandterm_onset = - Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, eventfields = [:onset]); -Xdc = Unfold.time_expand(X, timeexpandterm_onset, tbl) -@test_throws ArgumentError Unfold.time_expand(X, timeexpandterm_latency, tbl) + tbl2 = tbl = DataFrame([1 4]', [:onset]) + + timeexpandterm_latency = Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction) + timeexpandterm_onset = Unfold.TimeExpandedTerm( + FormulaTerm(Term, Term), + basisfunction, + eventfields = [:onset], + ) + Xdc = Unfold.time_expand(X, timeexpandterm_onset, tbl) + @test_throws ArgumentError Unfold.time_expand(X, timeexpandterm_latency, tbl) end @testset "combining designmatrices" begin -tbl = DataFrame([1 4]', [:latency]) -X = ones(size(tbl)) -basisfunction1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1") -basisfunction2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2") -f = @formula 0 ~ 1 -Xdc1 = designmatrix(UnfoldLinearModel, f, tbl, basisfunction1) -Xdc2 = designmatrix(UnfoldLinearModel, f, tbl .+ 1, basisfunction2) - -Xdc = Xdc1 + Xdc2 -@test size(modelmatrix(Xdc), 2) == size(modelmatrix(Xdc1), 2) + size(modelmatrix(Xdc2), 2) -@test length(Xdc.events) == 2 - -Xdc_3 = Xdc1 + Xdc2 + Xdc2 - -@test size(modelmatrix(Xdc_3), 2) == size(modelmatrix(Xdc1), 2) + 2*size(modelmatrix(Xdc2), 2) -@test length(Xdc_3.events) == 3 - - -basisfunction1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1") -basisfunction2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2") - -tbl = DataFrame( - [1 4 10 15 20 22 31 37; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]', - [:latency, :subject, :item], -) -tbl2 = DataFrame( - [2 3 12 18 19 25 40 43; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]', - [:latency, :subject, :itemB], -) -y = Float64.([collect(range(1, stop = 100))...])' -transform!(tbl, :subject => categorical => :subject) -transform!(tbl2, :itemB => categorical => :itemB) -transform!(tbl, :item => categorical => :item) -#tbl.itemB = tbl.item -f3 = @formula 0 ~ 1 + (1 | item) + (1 | subject) -f4 = @formula 0 ~ 1 + (1 | itemB) -f4_wrong = @formula 0 ~ 1 + (1 | item) -Xdc3 = designmatrix(UnfoldLinearMixedModel, f3, tbl, basisfunction1) -Xdc4 = designmatrix(UnfoldLinearMixedModel, f4, tbl2, basisfunction2) -Xdc4_wrong = designmatrix(UnfoldLinearMixedModel, f4_wrong, tbl, basisfunction2) - -Xdc = Xdc3 + Xdc4; -@test typeof(Xdc.Xs[1][1]) <: SparseArrays.SparseMatrixCSC -@test length(Xdc.Xs) == 4 # one FeMat + 3 ReMat -@test_throws String Xdc3 + Xdc4_wrong -uf = UnfoldLinearMixedModelContinuousTime(Dict(), Xdc, []) -Unfold.fit!(uf, y); - -Xdc = Xdc3 + Xdc4; -df = fit( - UnfoldLinearMixedModelContinuousTime, - Xdc, - rand(1, size(Unfold.modelmatrix(Xdc)[1], 1)), -) -@test size(Unfold.coef(df), 2) == 17 - -## Speedtest -if 1 == 0 - x = collect(range(10, stop = 10000, step = 10)) - tbl = DataFrame(reshape(x, 1000, 1), [:latency]) - X = ones(size(tbl, 1), 3) .* [1, 2, 3]' - - - basisfunction = firbasis(τ = (0, 1), sfreq = 100, name = "test") - term = Unfold.TimeExpandedTerm(Term, basisfunction, :latency) - @time Xdc = Matrix(Unfold.time_expand(X, term, tbl)) - -end + tbl = DataFrame([1 4]', [:latency]) + X = ones(size(tbl)) + basisfunction1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1") + basisfunction2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2") + f = @formula 0 ~ 1 + Xdc1 = designmatrix(UnfoldLinearModel, f, tbl, basisfunction1) + Xdc2 = designmatrix(UnfoldLinearModel, f, tbl .+ 1, basisfunction2) + + Xdc = Xdc1 + Xdc2 + @test size(modelmatrix(Xdc), 2) == + size(modelmatrix(Xdc1), 2) + size(modelmatrix(Xdc2), 2) + @test length(Xdc.events) == 2 + + Xdc_3 = Xdc1 + Xdc2 + Xdc2 + + @test size(modelmatrix(Xdc_3), 2) == + size(modelmatrix(Xdc1), 2) + 2 * size(modelmatrix(Xdc2), 2) + @test length(Xdc_3.events) == 3 + + + basisfunction1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1") + basisfunction2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2") + + tbl = DataFrame( + [1 4 10 15 20 22 31 37; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]', + [:latency, :subject, :item], + ) + tbl2 = DataFrame( + [2 3 12 18 19 25 40 43; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]', + [:latency, :subject, :itemB], + ) + y = Float64.([collect(range(1, stop = 100))...])' + transform!(tbl, :subject => categorical => :subject) + transform!(tbl2, :itemB => categorical => :itemB) + transform!(tbl, :item => categorical => :item) + #tbl.itemB = tbl.item + f3 = @formula 0 ~ 1 + (1 | item) + (1 | subject) + f4 = @formula 0 ~ 1 + (1 | itemB) + f4_wrong = @formula 0 ~ 1 + (1 | item) + Xdc3 = designmatrix(UnfoldLinearMixedModel, f3, tbl, basisfunction1) + Xdc4 = designmatrix(UnfoldLinearMixedModel, f4, tbl2, basisfunction2) + Xdc4_wrong = designmatrix(UnfoldLinearMixedModel, f4_wrong, tbl, basisfunction2) + + Xdc = Xdc3 + Xdc4 + @test typeof(Xdc.Xs[1][1]) <: SparseArrays.SparseMatrixCSC + @test length(Xdc.Xs) == 4 # one FeMat + 3 ReMat + @test_throws String Xdc3 + Xdc4_wrong + uf = UnfoldLinearMixedModelContinuousTime(Dict(), Xdc, []) + Unfold.fit!(uf, y) + + Xdc = Xdc3 + Xdc4 + df = fit( + UnfoldLinearMixedModelContinuousTime, + Xdc, + rand(1, size(Unfold.modelmatrix(Xdc)[1], 1)), + ) + @test size(Unfold.coef(df), 2) == 17 + + ## Speedtest + if 1 == 0 + x = collect(range(10, stop = 10000, step = 10)) + tbl = DataFrame(reshape(x, 1000, 1), [:latency]) + X = ones(size(tbl, 1), 3) .* [1, 2, 3]' + + + basisfunction = firbasis(τ = (0, 1), sfreq = 100, name = "test") + term = Unfold.TimeExpandedTerm(Term, basisfunction, :latency) + @time Xdc = Matrix(Unfold.time_expand(X, term, tbl)) + + end end @testset "equalizeReMatLengths" begin -bf1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1") -bf2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2") - -tbl1 = DataFrame( - [1 4 10 15 20 22 31 37; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]', - [:latency, :subject, :item], -) -tbl2 = DataFrame( - [2 3 12 18 19 25 40 43; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]', - [:latency, :subject, :itemB], -) - -transform!(tbl1, :subject => categorical => :subject) -transform!(tbl1, :item => categorical => :item) -transform!(tbl2, :itemB => categorical => :itemB) -#tbl.itemB = tbl.item -f1 = @formula 0 ~ 1 + (1 | item) + (1 | subject) -f2 = @formula 0 ~ 1 + (1 | itemB) - -form = apply_schema(f1, schema(f1, tbl1), MixedModels.LinearMixedModel) -form = Unfold.apply_basisfunction(form, bf1, nothing) -X1 = modelcols.(form.rhs, Ref(tbl1)) - -form = apply_schema(f2, schema(f2, tbl2), MixedModels.LinearMixedModel) -form = Unfold.apply_basisfunction(form, bf2, nothing) -X2 = modelcols.(form.rhs, Ref(tbl2)) - -# no missmatch, shouldnt change anything then -X = deepcopy(X1[2:end]) -if !isdefined(Base, :get_extension) - include("../ext/UnfoldMixedModelsExt/UnfoldMixedModelsExt.jl") - ext = UnfoldMixedModelsExt -else - ext = Base.get_extension(Unfold,:UnfoldMixedModelsExt) -end -ext.equalizeReMatLengths!(X) -@test all([x[1] for x in size.(X)] .== 48) - -X = (deepcopy(X1[2:end])..., deepcopy(X2[2:end])...) -@test !all([x[1] for x in size.(X)] .== 49) # not alllenghts the same -ext.equalizeReMatLengths!(X) -@test all([x[1] for x in size.(X)] .== 49) # now all lengths the same :-) - - -X = deepcopy(X2[2]) - -@test size(X)[1] == 49 -ext.changeReMatSize!(X, 52) -@test size(X)[1] == 52 - -X = deepcopy(X2[2]) -@test size(X)[1] == 49 -ext.changeReMatSize!(X, 40) -@test size(X)[1] == 40 - - -X = (deepcopy(X1)..., deepcopy(X2[2:end])...) -@test size(X[1])[1] == 48 -@test size(X[2])[1] == 48 -@test size(X[3])[1] == 48 -@test size(X[4])[1] == 49 -XA, XB = ext.changeMatSize!(52, X[1], X[2:end]) -@test size(XA)[1] == 52 -@test size(XB)[1] == 52 - -XA, XB = ext.changeMatSize!(40, X[1], X[2:end]) -@test size(XA)[1] == 40 -@test size(XB)[1] == 40 - -XA, XB = ext.changeMatSize!(30, Matrix(X[1]), X[2:end]) -@test size(XA)[1] == 30 -@test size(XB)[1] == 30 + bf1 = firbasis(τ = (0, 1), sfreq = 10, name = "basis1") + bf2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2") + + tbl1 = DataFrame( + [1 4 10 15 20 22 31 37; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]', + [:latency, :subject, :item], + ) + tbl2 = DataFrame( + [2 3 12 18 19 25 40 43; 1 1 1 2 2 2 3 3; 1 2 3 1 2 3 1 2]', + [:latency, :subject, :itemB], + ) + + transform!(tbl1, :subject => categorical => :subject) + transform!(tbl1, :item => categorical => :item) + transform!(tbl2, :itemB => categorical => :itemB) + #tbl.itemB = tbl.item + f1 = @formula 0 ~ 1 + (1 | item) + (1 | subject) + f2 = @formula 0 ~ 1 + (1 | itemB) + + form = apply_schema(f1, schema(f1, tbl1), MixedModels.LinearMixedModel) + form = Unfold.apply_basisfunction(form, bf1, nothing) + X1 = modelcols.(form.rhs, Ref(tbl1)) + + form = apply_schema(f2, schema(f2, tbl2), MixedModels.LinearMixedModel) + form = Unfold.apply_basisfunction(form, bf2, nothing) + X2 = modelcols.(form.rhs, Ref(tbl2)) + + # no missmatch, shouldnt change anything then + X = deepcopy(X1[2:end]) + if !isdefined(Base, :get_extension) + include("../ext/UnfoldMixedModelsExt/UnfoldMixedModelsExt.jl") + ext = UnfoldMixedModelsExt + else + ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) + end + ext.equalizeReMatLengths!(X) + @test all([x[1] for x in size.(X)] .== 48) + + X = (deepcopy(X1[2:end])..., deepcopy(X2[2:end])...) + @test !all([x[1] for x in size.(X)] .== 49) # not alllenghts the same + ext.equalizeReMatLengths!(X) + @test all([x[1] for x in size.(X)] .== 49) # now all lengths the same :-) + + + X = deepcopy(X2[2]) + + @test size(X)[1] == 49 + ext.changeReMatSize!(X, 52) + @test size(X)[1] == 52 + + X = deepcopy(X2[2]) + @test size(X)[1] == 49 + ext.changeReMatSize!(X, 40) + @test size(X)[1] == 40 + + + X = (deepcopy(X1)..., deepcopy(X2[2:end])...) + @test size(X[1])[1] == 48 + @test size(X[2])[1] == 48 + @test size(X[3])[1] == 48 + @test size(X[4])[1] == 49 + XA, XB = ext.changeMatSize!(52, X[1], X[2:end]) + @test size(XA)[1] == 52 + @test size(XB)[1] == 52 + + XA, XB = ext.changeMatSize!(40, X[1], X[2:end]) + @test size(XA)[1] == 40 + @test size(XB)[1] == 40 + + XA, XB = ext.changeMatSize!(30, Matrix(X[1]), X[2:end]) + @test size(XA)[1] == 30 + @test size(XB)[1] == 30 end @testset "Some LinearMixedModel tests" begin -data, evts = loadtestdata("testCase3", dataPath = (@__DIR__) * "/data") # -evts.subject = categorical(evts.subject) + data, evts = loadtestdata("testCase3", dataPath = (@__DIR__) * "/data") # + evts.subject = categorical(evts.subject) -f_zc = @formula 0 ~ 1 + condA + condB + zerocorr(1 + condA + condB | subject) -basisfunction = firbasis(τ = (-0.1, 0.1), sfreq = 10, name = "ABC") -Xdc_zc = designmatrix(UnfoldLinearMixedModel, f_zc, evts, basisfunction) + f_zc = @formula 0 ~ 1 + condA + condB + zerocorr(1 + condA + condB | subject) + basisfunction = firbasis(τ = (-0.1, 0.1), sfreq = 10, name = "ABC") + Xdc_zc = designmatrix(UnfoldLinearMixedModel, f_zc, evts, basisfunction) -@test length(Xdc_zc.Xs[2].inds) == 9 -f = @formula 0 ~ 1 + condA + condB + (1 + condA + condB | subject) -Xdc = designmatrix(UnfoldLinearMixedModel, f, evts, basisfunction) -@test length(Xdc.Xs[2].inds) == (9 * 9 + 9) / 2 + @test length(Xdc_zc.Xs[2].inds) == 9 + f = @formula 0 ~ 1 + condA + condB + (1 + condA + condB | subject) + Xdc = designmatrix(UnfoldLinearMixedModel, f, evts, basisfunction) + @test length(Xdc.Xs[2].inds) == (9 * 9 + 9) / 2 -# test bug with not sequential subjects -evts_nonseq = copy(evts) -evts_nonseq = evts_nonseq[.!(evts_nonseq.subject .== 2), :] -Xdc_nonseq = designmatrix(UnfoldLinearMixedModel, f_zc, evts_nonseq, basisfunction) -# This used to lead to problems here: -fit(UnfoldLinearMixedModel, Xdc_nonseq, data'); + # test bug with not sequential subjects + evts_nonseq = copy(evts) + evts_nonseq = evts_nonseq[.!(evts_nonseq.subject .== 2), :] + Xdc_nonseq = designmatrix(UnfoldLinearMixedModel, f_zc, evts_nonseq, basisfunction) + # This used to lead to problems here: + fit(UnfoldLinearMixedModel, Xdc_nonseq, data') end #basisfunction2 = firbasis(τ = (0, 0.5), sfreq = 10, name = "basis2") @testset "Missings in Events" begin -tbl = DataFrame( - :a => [1,2,3,4,5,6,7,8], - :b=>[1,1,1,2,2,2,3,missing], - :c=>[1,2,3,4,5,6,7,missing], - :d=>["1","2","3","4","5","6","7","8"], - :e=>["1","2","3","4","5","6","7",missing], - :event=>[1,1,1,1,2,2,2,2], - :latency=> [10,20,30,40,50,60,70,80]) -tbl.event = string.(tbl.event) - designmatrix(UnfoldLinearModel,@formula(0~a),tbl) - @test_throws ErrorException designmatrix(UnfoldLinearModel,@formula(0~a+b),tbl) - @test_throws ErrorException designmatrix(UnfoldLinearModel,@formula(0~e),tbl) + tbl = DataFrame( + :a => [1, 2, 3, 4, 5, 6, 7, 8], + :b => [1, 1, 1, 2, 2, 2, 3, missing], + :c => [1, 2, 3, 4, 5, 6, 7, missing], + :d => ["1", "2", "3", "4", "5", "6", "7", "8"], + :e => ["1", "2", "3", "4", "5", "6", "7", missing], + :event => [1, 1, 1, 1, 2, 2, 2, 2], + :latency => [10, 20, 30, 40, 50, 60, 70, 80], + ) + tbl.event = string.(tbl.event) + designmatrix(UnfoldLinearModel, @formula(0 ~ a), tbl) + @test_throws ErrorException designmatrix(UnfoldLinearModel, @formula(0 ~ a + b), tbl) + @test_throws ErrorException designmatrix(UnfoldLinearModel, @formula(0 ~ e), tbl) # including an actual missing doesnt work - design = Dict("1"=>(@formula(0~a+b+c+d+e),firbasis((0,1),1)),"2"=>(@formula(0~a+b+c+d+e),firbasis((0,1),1))) + design = Dict( + "1" => (@formula(0 ~ a + b + c + d + e), firbasis((0, 1), 1)), + "2" => (@formula(0 ~ a + b + c + d + e), firbasis((0, 1), 1)), + ) uf = UnfoldLinearModelContinuousTime(design) - @test_throws ErrorException designmatrix(uf,tbl); + @test_throws ErrorException designmatrix(uf, tbl) # but if the missing is in another event, no problem - design = Dict("1"=>(@formula(0~a+b+c+d+e),firbasis((0,1),1)),"2"=>(@formula(0~a+d),firbasis((0,1),1))) + design = Dict( + "1" => (@formula(0 ~ a + b + c + d + e), firbasis((0, 1), 1)), + "2" => (@formula(0 ~ a + d), firbasis((0, 1), 1)), + ) uf = UnfoldLinearModelContinuousTime(design) - designmatrix(uf,tbl); + designmatrix(uf, tbl) # prior to the Missing disallow sanity check, this gave an error - design = Dict("1"=>(@formula(0~spl(a,4)+spl(b,4)+d+e),firbasis((0,1),1)),"2"=>(@formula(0~a+d),firbasis((0,1),1))) + design = Dict( + "1" => (@formula(0 ~ spl(a, 4) + spl(b, 4) + d + e), firbasis((0, 1), 1)), + "2" => (@formula(0 ~ a + d), firbasis((0, 1), 1)), + ) uf = UnfoldLinearModelContinuousTime(design) - designmatrix(uf,tbl); -end \ No newline at end of file + designmatrix(uf, tbl) +end diff --git a/test/effects.jl b/test/effects.jl index d67f735f..457d55cd 100644 --- a/test/effects.jl +++ b/test/effects.jl @@ -4,28 +4,31 @@ data_r = reshape(data, (1, :)) data_e, times = Unfold.epoch(data = data_r, tbl = evts, τ = (0, 0.05), sfreq = 10) # cut the data into epochs f = @formula 0 ~ 1 + conditionA + continuousA # 1 -m_mul = fit(Unfold.UnfoldModel, Dict(Any=>(f,times)), evts, data_e) +m_mul = fit(Unfold.UnfoldModel, Dict(Any => (f, times)), evts, data_e) ## @testset "Mass Univariate, all specified" begin - # test simple case - eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),m_mul) - @test size(eff,1) == 2 # we specified 2 levels @ 1 time point - @test eff.conditionA ≈ [0.,1.] # we want to different levels - @test eff.yhat ≈ [2.0,5.0] # these are the perfect predicted values + # test simple case + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [0]), m_mul) + @test size(eff, 1) == 2 # we specified 2 levels @ 1 time point + @test eff.conditionA ≈ [0.0, 1.0] # we want to different levels + @test eff.yhat ≈ [2.0, 5.0] # these are the perfect predicted values - # combination 2 levels / 6 values - eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[-2,0,2]),m_mul) - @test size(eff,1) == 6 # we want 6 values - @test eff.conditionA ≈ [0.,0.,0.,1.,1.,1.] - @test eff.continuousA ≈ [-2,0,2,-2,0,2.] + # combination 2 levels / 6 values + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [-2, 0, 2]), m_mul) + @test size(eff, 1) == 6 # we want 6 values + @test eff.conditionA ≈ [0.0, 0.0, 0.0, 1.0, 1.0, 1.0] + @test eff.continuousA ≈ [-2, 0, 2, -2, 0, 2.0] end @testset "Mass Univariate, typified" begin -# testing typical value - eff_man = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[mean(evts.continuousA)]),m_mul) - eff_typ = Unfold.effects(Dict(:conditionA => [0,1]),m_mul) - @test eff_man.yhat ≈ eff_typ.yhat + # testing typical value + eff_man = Unfold.effects( + Dict(:conditionA => [0, 1], :continuousA => [mean(evts.continuousA)]), + m_mul, + ) + eff_typ = Unfold.effects(Dict(:conditionA => [0, 1]), m_mul) + @test eff_man.yhat ≈ eff_typ.yhat end ## Testing Splines @@ -34,21 +37,24 @@ m_mul_spl = fit(UnfoldModel, f_spl, evts, data_e, times) @testset "Mass Univariate, splines" begin -eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),m_mul_spl) -@test size(eff,1) == 2 # we specified 2 levels @ 1 time point -@test eff.conditionA ≈ [0.,1.] # we want to different levels -@test eff.yhat ≈ [2.0,5.0] # these are the perfect predicted values - -# combination 2 levels / 6 values -eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[-0.5,0,0.5]),m_mul_spl) -@test size(eff,1) == 6 # we want 6 values -@test eff.conditionA ≈ [0.,0.,0.,1.,1.,1.] -@test eff.continuousA ≈ [-0.5,0,0.5,-0.5,0,0.5] -@test eff.yhat ≈ [0,2,4,3,5,7] - -# testing for safe predictions -eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[2]),m_mul_spl) -@test all(ismissing.(eff.yhat )) + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [0]), m_mul_spl) + @test size(eff, 1) == 2 # we specified 2 levels @ 1 time point + @test eff.conditionA ≈ [0.0, 1.0] # we want to different levels + @test eff.yhat ≈ [2.0, 5.0] # these are the perfect predicted values + + # combination 2 levels / 6 values + eff = Unfold.effects( + Dict(:conditionA => [0, 1], :continuousA => [-0.5, 0, 0.5]), + m_mul_spl, + ) + @test size(eff, 1) == 6 # we want 6 values + @test eff.conditionA ≈ [0.0, 0.0, 0.0, 1.0, 1.0, 1.0] + @test eff.continuousA ≈ [-0.5, 0, 0.5, -0.5, 0, 0.5] + @test eff.yhat ≈ [0, 2, 4, 3, 5, 7] + + # testing for safe predictions + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [2]), m_mul_spl) + @test all(ismissing.(eff.yhat)) end ## Timeexpansion @@ -57,12 +63,12 @@ f = @formula 0 ~ 1 + conditionA + continuousA # 1 @testset "Time Expansion, one event" begin -uf = fit(Unfold.UnfoldModel, Dict(Any=>(f,firbasis([0,0.1],10))), evts, data) -eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),uf) -@test nrow(eff) == 4 -@test eff.yhat ≈ [2., 2., 5., 5.] -@test eff.conditionA ≈ [0.,0.,1.,1.] -@test eff.continuousA ≈ [0.,0.,0.,0.] + uf = fit(Unfold.UnfoldModel, Dict(Any => (f, firbasis([0, 0.1], 10))), evts, data) + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [0]), uf) + @test nrow(eff) == 4 + @test eff.yhat ≈ [2.0, 2.0, 5.0, 5.0] + @test eff.conditionA ≈ [0.0, 0.0, 1.0, 1.0] + @test eff.continuousA ≈ [0.0, 0.0, 0.0, 0.0] end data, evts = loadtestdata("test_case_4a") # @@ -70,149 +76,191 @@ b1 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisA") b2 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisB") b3 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisC") f = @formula 0 ~ 1 # 1 -m_tul = fit(UnfoldModel, Dict("eventA"=>(f,b1),"eventB"=>(f,b2)), evts, data,eventcolumn="type") +m_tul = fit( + UnfoldModel, + Dict("eventA" => (f, b1), "eventB" => (f, b2)), + evts, + data, + eventcolumn = "type", +) @testset "Time Expansion, two events" begin -eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),m_tul) -@test unique(eff.basisname)==["basisA","basisB"] -@test unique(eff.yhat) ≈ [2,3] -@test size(eff,1) == 2*2*20 # 2 basisfunctions, 2x conditionA, 1s a 20Hz + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [0]), m_tul) + @test unique(eff.basisname) == ["basisA", "basisB"] + @test unique(eff.yhat) ≈ [2, 3] + @test size(eff, 1) == 2 * 2 * 20 # 2 basisfunctions, 2x conditionA, 1s a 20Hz -eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[-1,0,1]),m_tul) -@test size(eff,1) == 2*6*20 + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [-1, 0, 1]), m_tul) + @test size(eff, 1) == 2 * 6 * 20 end @testset "Time Expansion, three events" begin evts_3 = deepcopy(evts) - evts_3.type[1:50] .= "eventC" - evts_3.type[50:100] .= "eventD" - m_tul_3 = fit(UnfoldModel, Dict("eventA"=>(f,b1),"eventB"=>(f,b2),"eventC"=>(f,b3)), evts_3, data,eventcolumn="type") - - eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[0]),m_tul_3) - - @test size(eff,1) == 3*2*20 # 2 basisfunctions, 2x conditionA, 1s a 20Hz - + evts_3.type[1:50] .= "eventC" + evts_3.type[50:100] .= "eventD" + m_tul_3 = fit( + UnfoldModel, + Dict("eventA" => (f, b1), "eventB" => (f, b2), "eventC" => (f, b3)), + evts_3, + data, + eventcolumn = "type", + ) + + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [0]), m_tul_3) + + @test size(eff, 1) == 3 * 2 * 20 # 2 basisfunctions, 2x conditionA, 1s a 20Hz + end @testset "Time Expansion, two events different size + different formulas" begin -## Different sized events + different Formulas -data, evts = loadtestdata("test_case_4a") # -evts[!,:continuousA] = rand(MersenneTwister(42),nrow(evts)) -b1 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisA") -b2 = firbasis(τ = (0.0, 0.5), sfreq = 20, name = "basisB") -f1 = @formula 0 ~ 1 # 1 -f2 = @formula 0 ~ 1+continuousA # 1 -m_tul = fit(UnfoldModel, Dict("eventA"=>(f1,b1),"eventB"=>(f2,b2)), evts, data,eventcolumn="type") -eff = Unfold.effects(Dict(:conditionA => [0,1],:continuousA =>[-1,0,1]),m_tul) -@test nrow(eff) == (length(Unfold.times(b1))-1+length(Unfold.times(b2))-1)*6 -@test sum(eff.basisname .== "basisA") == 120 -@test sum(eff.basisname .== "basisB") == 66 + ## Different sized events + different Formulas + data, evts = loadtestdata("test_case_4a") # + evts[!, :continuousA] = rand(MersenneTwister(42), nrow(evts)) + b1 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisA") + b2 = firbasis(τ = (0.0, 0.5), sfreq = 20, name = "basisB") + f1 = @formula 0 ~ 1 # 1 + f2 = @formula 0 ~ 1 + continuousA # 1 + m_tul = fit( + UnfoldModel, + Dict("eventA" => (f1, b1), "eventB" => (f2, b2)), + evts, + data, + eventcolumn = "type", + ) + eff = Unfold.effects(Dict(:conditionA => [0, 1], :continuousA => [-1, 0, 1]), m_tul) + @test nrow(eff) == (length(Unfold.times(b1)) - 1 + length(Unfold.times(b2)) - 1) * 6 + @test sum(eff.basisname .== "basisA") == 120 + @test sum(eff.basisname .== "basisB") == 66 end ## Test two channels data, evts = loadtestdata("test_case_3a") # -data_r = repeat(reshape(data, (1, :)),3,1) -data_r[2,:] = data_r[2,:] .* 2 +data_r = repeat(reshape(data, (1, :)), 3, 1) +data_r[2, :] = data_r[2, :] .* 2 data_e, times = Unfold.epoch(data = data_r, tbl = evts, τ = (0, 0.05), sfreq = 10) # cut the data into epochs # f = @formula 0 ~ 1 + conditionA + continuousA # 1 -m_mul = fit(Unfold.UnfoldModel, Dict(Any=>(f,times)), evts, data_e) -m_tul = fit(Unfold.UnfoldModel, Dict(Any=>(f,firbasis([0,.05],10))), evts, data_r) +m_mul = fit(Unfold.UnfoldModel, Dict(Any => (f, times)), evts, data_e) +m_tul = fit(Unfold.UnfoldModel, Dict(Any => (f, firbasis([0, 0.05], 10))), evts, data_r) @testset "Two channels" begin - # test simple case - eff_m = Unfold.effects(Dict(:conditionA => [0,1,0,1],:continuousA =>[0]),m_mul) - eff_t = Unfold.effects(Dict(:conditionA => [0,1,0,1],:continuousA =>[0]),m_tul) + # test simple case + eff_m = Unfold.effects(Dict(:conditionA => [0, 1, 0, 1], :continuousA => [0]), m_mul) + eff_t = Unfold.effects(Dict(:conditionA => [0, 1, 0, 1], :continuousA => [0]), m_tul) - @test eff_m.yhat ≈ eff_t.yhat - @test length(unique(eff_m.channel)) == 3 - @test eff_m[eff_m.channel .==1,:yhat] ≈ eff_m[eff_m.channel .==2,:yhat]./2 - @test eff_m[eff_m.channel .==1,:yhat] ≈ [2,5,2,5.] # these are the perfect predicted values - note that we requested them twice + @test eff_m.yhat ≈ eff_t.yhat + @test length(unique(eff_m.channel)) == 3 + @test eff_m[eff_m.channel.==1, :yhat] ≈ eff_m[eff_m.channel.==2, :yhat] ./ 2 + @test eff_m[eff_m.channel.==1, :yhat] ≈ [2, 5, 2, 5.0] # these are the perfect predicted values - note that we requested them twice end @testset "Timeexpansion, two events, typified" begin -data, evts = loadtestdata("test_case_4a") # -evts[!,:continuousA] = rand(MersenneTwister(42),nrow(evts)) -evts[!,:continuousB] = rand(MersenneTwister(43),nrow(evts)) -ixA = evts.type .== "eventA" -evts.continuousB[ixA] = evts.continuousB[ixA] .-mean(evts.continuousB[ixA]) .-5 -evts.continuousB[.!ixA] = evts.continuousB[.!ixA] .-mean(evts.continuousB[.!ixA]) .+ 0.5 -b1 = firbasis(τ = (0.0, 0.02), sfreq = 20, name = "basisA") -b2 = firbasis(τ = (1.0, 1.02), sfreq = 20, name = "basisB") -f1 = @formula 0 ~ 1+continuousA # 1 -f2 = @formula 0 ~ 1+continuousB # 1 -m_tul = fit(UnfoldModel, Dict("eventA"=>(f1,b1),"eventB"=>(f2,b2)), evts, data,eventcolumn="type") - -m_tul.modelfit.estimate .= [0 -1 0 6] -eff = Unfold.effects(Dict(:continuousA => [0,1]),m_tul) -eff = Unfold.effects(Dict(:continuousA => [0,1],:continuousB =>[0.5]),m_tul) - - -@test eff.yhat[3] == eff.yhat[4] -@test eff.yhat[1] == 0. -@test eff.yhat[2] == -1. -@test eff.yhat[3] == 3 + data, evts = loadtestdata("test_case_4a") # + evts[!, :continuousA] = rand(MersenneTwister(42), nrow(evts)) + evts[!, :continuousB] = rand(MersenneTwister(43), nrow(evts)) + ixA = evts.type .== "eventA" + evts.continuousB[ixA] = evts.continuousB[ixA] .- mean(evts.continuousB[ixA]) .- 5 + evts.continuousB[.!ixA] = + evts.continuousB[.!ixA] .- mean(evts.continuousB[.!ixA]) .+ 0.5 + b1 = firbasis(τ = (0.0, 0.02), sfreq = 20, name = "basisA") + b2 = firbasis(τ = (1.0, 1.02), sfreq = 20, name = "basisB") + f1 = @formula 0 ~ 1 + continuousA # 1 + f2 = @formula 0 ~ 1 + continuousB # 1 + m_tul = fit( + UnfoldModel, + Dict("eventA" => (f1, b1), "eventB" => (f2, b2)), + evts, + data, + eventcolumn = "type", + ) + + m_tul.modelfit.estimate .= [0 -1 0 6] + eff = Unfold.effects(Dict(:continuousA => [0, 1]), m_tul) + eff = Unfold.effects(Dict(:continuousA => [0, 1], :continuousB => [0.5]), m_tul) + + + @test eff.yhat[3] == eff.yhat[4] + @test eff.yhat[1] == 0.0 + @test eff.yhat[2] == -1.0 + @test eff.yhat[3] == 3 end @testset "timeexpansion, Interactions, two events" begin - data, evts = loadtestdata("test_case_4a") # - evts[!,:continuousA] = rand(MersenneTwister(42),nrow(evts)) - evts[!,:continuousB] = ["m","x"][Int.(1 .+ round.(rand(MersenneTwister(43),nrow(evts))))] - + data, evts = loadtestdata("test_case_4a") # + evts[!, :continuousA] = rand(MersenneTwister(42), nrow(evts)) + evts[!, :continuousB] = + ["m", "x"][Int.(1 .+ round.(rand(MersenneTwister(43), nrow(evts))))] - b1 = firbasis(τ = (0.0, 0.02), sfreq = 20, name = "basisA") - b2 = firbasis(τ = (1.0, 1.02), sfreq = 20, name = "basisB") - f1 = @formula 0 ~ 1+continuousA*continuousB # 1 - f2 = @formula 0 ~ 1+continuousB # 1 - m_tul = fit(UnfoldModel, Dict("eventA"=>(f1,b1),"eventB"=>(f2,b2)), evts, data,eventcolumn="type") - m_tul.modelfit.estimate .= [0, -1, 0, 2., 0.,0.]' + b1 = firbasis(τ = (0.0, 0.02), sfreq = 20, name = "basisA") + b2 = firbasis(τ = (1.0, 1.02), sfreq = 20, name = "basisB") + f1 = @formula 0 ~ 1 + continuousA * continuousB # 1 + f2 = @formula 0 ~ 1 + continuousB # 1 + m_tul = fit( + UnfoldModel, + Dict("eventA" => (f1, b1), "eventB" => (f2, b2)), + evts, + data, + eventcolumn = "type", + ) + + m_tul.modelfit.estimate .= [0, -1, 0, 2.0, 0.0, 0.0]' - eff = Unfold.effects(Dict(:continuousA => [0,1]),m_tul) - @test size(eff,1) == 4 - @test all(eff.basisname[1:2] .== "basisA") - @test all(eff.basisname[4:end] .== "basisB") - @test eff.yhat ≈ [0., mean(evts.continuousB[evts.type .== "eventA"].=="x") * coef(m_tul)[4] + 1*coef(m_tul)[2], 0., 0.] + eff = Unfold.effects(Dict(:continuousA => [0, 1]), m_tul) + @test size(eff, 1) == 4 + @test all(eff.basisname[1:2] .== "basisA") + @test all(eff.basisname[4:end] .== "basisB") + @test eff.yhat ≈ [ + 0.0, + mean(evts.continuousB[evts.type.=="eventA"] .== "x") * coef(m_tul)[4] + + 1 * coef(m_tul)[2], + 0.0, + 0.0, + ] - eff = Unfold.effects(Dict(:continuousB => ["m","x"]),m_tul) - @test eff.yhat[1] == -eff.yhat[2] - @test all(eff.yhat[3:4] .≈ 0.) + eff = Unfold.effects(Dict(:continuousB => ["m", "x"]), m_tul) + @test eff.yhat[1] == -eff.yhat[2] + @test all(eff.yhat[3:4] .≈ 0.0) - eff = Unfold.effects(Dict(:continuousA => [0,1],:continuousB => ["m","x"]),m_tul) - @test eff.yhat[3:4] == [-1,1] - @test all(eff.yhat[1:2, 5:end] .== 0) + eff = Unfold.effects(Dict(:continuousA => [0, 1], :continuousB => ["m", "x"]), m_tul) + @test eff.yhat[3:4] == [-1, 1] + @test all(eff.yhat[1:2, 5:end] .== 0) end @testset "splinesMissings" begin - @testset "Missings in Events" begin - tbl = DataFrame( - :a => [1,2,3,4,5,6,7,8], - :b=>[1,1,1,2,2,2,3,missing], - :c=>[1,2,3,4,5,6,7,missing], - :d=>["1","2","3","4","5","6","7","8"], - :e=>["1","2","3","4","5","6","7",missing], - :event=>[1,1,1,1,2,2,2,2], - :latency=> [10,20,30,40,50,60,70,80]) - data = rand(120) - tbl.event = string.(tbl.event) - - # prior to the Missing disallow sanity check, this gave an error - design = Dict("1"=>(@formula(0~spl(a,4)+spl(b,4)+d+e),firbasis((0,2),1)),"2"=>(@formula(0~a+d),firbasis((0,1),1))) - m =fit(UnfoldModel,design,tbl,data) - effects(Dict(:a=>[1,2,3],:b=>[1,1.5]),m) - end + @testset "Missings in Events" begin + tbl = DataFrame( + :a => [1, 2, 3, 4, 5, 6, 7, 8], + :b => [1, 1, 1, 2, 2, 2, 3, missing], + :c => [1, 2, 3, 4, 5, 6, 7, missing], + :d => ["1", "2", "3", "4", "5", "6", "7", "8"], + :e => ["1", "2", "3", "4", "5", "6", "7", missing], + :event => [1, 1, 1, 1, 2, 2, 2, 2], + :latency => [10, 20, 30, 40, 50, 60, 70, 80], + ) + data = rand(120) + tbl.event = string.(tbl.event) + + # prior to the Missing disallow sanity check, this gave an error + design = Dict( + "1" => (@formula(0 ~ spl(a, 4) + spl(b, 4) + d + e), firbasis((0, 2), 1)), + "2" => (@formula(0 ~ a + d), firbasis((0, 1), 1)), + ) + m = fit(UnfoldModel, design, tbl, data) + effects(Dict(:a => [1, 2, 3], :b => [1, 1.5]), m) + end end @@ -244,4 +292,4 @@ end data, ) @test_broken eff = effects(Dict(:condition => ["car", "face"]), m) -end \ No newline at end of file +end diff --git a/test/fit.jl b/test/fit.jl index a46e7e9d..3d67cc39 100644 --- a/test/fit.jl +++ b/test/fit.jl @@ -22,33 +22,43 @@ data_e, times = Unfold.epoch(data = data_r, tbl = evts, τ = (-1.0, 1.9), sfreq end @testset "epoched auto multi-event" begin - - evts_local = deepcopy(evts) - evts_local.type .= repeat(["A","B"],nrow(evts)÷2) - uf = fit(UnfoldModel, Dict("A" => (f, times)),evts_local,data_e;eventcolumn="type") - @test size(coef(uf)) == (2,59,3) - uf_2events = fit(UnfoldModel, Dict("A" => (f, times),"B"=>(@formula(0~1),times)),evts_local,data_e;eventcolumn="type") - @test size(coef(uf_2events)) == (2,59,4) + evts_local = deepcopy(evts) + evts_local.type .= repeat(["A", "B"], nrow(evts) ÷ 2) + + uf = fit(UnfoldModel, Dict("A" => (f, times)), evts_local, data_e; eventcolumn = "type") + @test size(coef(uf)) == (2, 59, 3) + uf_2events = fit( + UnfoldModel, + Dict("A" => (f, times), "B" => (@formula(0 ~ 1), times)), + evts_local, + data_e; + eventcolumn = "type", + ) + @test size(coef(uf_2events)) == (2, 59, 4) c = coeftable(uf) c2 = coeftable(uf_2events) - @test c2[c2.basisname .== "event: A",:] == c + @test c2[c2.basisname.=="event: A", :] == c + + e_uf2 = effects(Dict(:condtionA => [0, 1]), uf_2events) + e_uf = effects(Dict(:condtionA => [0, 1]), uf) - e_uf2 = effects(Dict(:condtionA=>[0,1]),uf_2events) - e_uf = effects(Dict(:condtionA=>[0,1]),uf) - end @testset "test Autodetection" begin - @test Unfold.designToModeltype(Dict(Any => (@formula(0 ~ 1), 0:10))) == UnfoldLinearModel + @test Unfold.designToModeltype(Dict(Any => (@formula(0 ~ 1), 0:10))) == + UnfoldLinearModel @test Unfold.designToModeltype(Dict(Any => (@formula(0 ~ 1 + A), 0:10))) == - UnfoldLinearModel + UnfoldLinearModel @test Unfold.designToModeltype( - Dict(Any => (@formula(0 ~ 1 + A), firbasis(τ = (-1, 1), sfreq = 20, name = "basisA"))), + Dict( + Any => + (@formula(0 ~ 1 + A), firbasis(τ = (-1, 1), sfreq = 20, name = "basisA")), + ), ) == UnfoldLinearModelContinuousTime @test Unfold.designToModeltype(Dict(Any => (@formula(0 ~ 1 + (1 | test)), 0:10))) == - UnfoldLinearMixedModel + UnfoldLinearMixedModel @test Unfold.designToModeltype( Dict( Any => ( @@ -62,9 +72,24 @@ end @testset "Bad Input" begin # check that if UnfoldLinearModel or UnfoldLinearModelContinuousTime is defined, that the design is appropriate basisfunction = firbasis(τ = (-1, 1), sfreq = 20, name = "basisA") - @test_throws "InputError" fit(UnfoldLinearModel,Dict(Any => (@formula(0 ~ 1),basisfunction)),evts,data_r) - @test_throws "InputError" fit(UnfoldLinearModel,Dict(Any => (@formula(0 ~ 1),basisfunction)),evts,data_e) - @test_throws "InputError" fit(UnfoldLinearModelContinuousTime,Dict(Any => (@formula(0 ~ 1),0:0.1:1)),evts,data_r) + @test_throws "InputError" fit( + UnfoldLinearModel, + Dict(Any => (@formula(0 ~ 1), basisfunction)), + evts, + data_r, + ) + @test_throws "InputError" fit( + UnfoldLinearModel, + Dict(Any => (@formula(0 ~ 1), basisfunction)), + evts, + data_e, + ) + @test_throws "InputError" fit( + UnfoldLinearModelContinuousTime, + Dict(Any => (@formula(0 ~ 1), 0:0.1:1)), + evts, + data_r, + ) end end @@ -72,12 +97,13 @@ end @testset "automatic, non-dictionary call" begin times = -1.0:0.05:1.9 m_mul = coeftable(fit(UnfoldLinearModel, f, evts, data_e, times)) - + @test m_mul[(m_mul.channel.==1).&(m_mul.time.==0.1), :estimate] ≈ [2, 3, 4] - - data_e_noreshape, times = Unfold.epoch(data = data, tbl = evts, τ = (-1.0, 1.9), sfreq = 20) # cut the data into epochs + + data_e_noreshape, times = + Unfold.epoch(data = data, tbl = evts, τ = (-1.0, 1.9), sfreq = 20) # cut the data into epochs m_mul_noreshape = coeftable(fit(UnfoldLinearModel, f, evts, data_e_noreshape, times)) @test m_mul_noreshape[ @@ -86,7 +112,8 @@ end ] ≈ [2, 3, 4] @test size(m_mul_noreshape)[1] == size(m_mul)[1] / 2 # test 2D call for UnfoldLinearModel - m_mul_autoreshape = coeftable(fit(UnfoldLinearModel, f, evts, data_e_noreshape[1,:,:], times)) + m_mul_autoreshape = + coeftable(fit(UnfoldLinearModel, f, evts, data_e_noreshape[1, :, :], times)) m_mul_autoreshape == m_mul_noreshape # Add Missing in Data data_e_missing = deepcopy(data_e) @@ -97,15 +124,18 @@ end # Special solver solver_lsmr_se with Standard Error se_solver = solver = (x, y) -> Unfold.solver_default(x, y, stderror = true) - m_mul_se = coeftable(Unfold.fit(UnfoldModel, f, evts, data_e, times, solver = se_solver)) + m_mul_se = + coeftable(Unfold.fit(UnfoldModel, f, evts, data_e, times, solver = se_solver)) @test all(m_mul_se.estimate .≈ m_mul.estimate) @test !all(isnothing.(m_mul_se.stderror)) # robust solver rob_solver = (x, y) -> Unfold.solver_robust(x, y)#,rlmOptions=(initial_coef=zeros(3 *length(times)),)) data_outlier = copy(data_e) - data_outlier[:,31,1] .= 1000 - m_mul_rob = coeftable(Unfold.fit(UnfoldModel, f, evts, data_outlier, times, solver = rob_solver)) + data_outlier[:, 31, 1] .= 1000 + m_mul_rob = coeftable( + Unfold.fit(UnfoldModel, f, evts, data_outlier, times, solver = rob_solver), + ) ix = findall(m_mul_rob.time .≈ 0.5) @test all(m_mul_rob.estimate[ix] .≈ m_mul.estimate[ix]) @@ -118,7 +148,7 @@ end basisfunction = firbasis(τ = (-1, 1), sfreq = 20, name = "basisA") @testset "timeexpanded univariate linear+missings" begin - + m_tul = coeftable(fit(UnfoldModel, f, evts, data_r, basisfunction)) @test isapprox( @@ -128,7 +158,7 @@ basisfunction = firbasis(τ = (-1, 1), sfreq = 20, name = "basisA") ) # test without reshape, i.e. 1 channel vector e.g. size(data) = (1200,) - m_tul_noreshape = coeftable(fit(UnfoldModel, f, evts, data, basisfunction)); + m_tul_noreshape = coeftable(fit(UnfoldModel, f, evts, data, basisfunction)) @test size(m_tul_noreshape)[1] == size(m_tul)[1] / 2 # Test under missing data @@ -182,7 +212,7 @@ basisfunction = firbasis(τ = (-1, 1), sfreq = 20, name = "basisA") end @testset "runtime tests" begin - + # runntime tests - does something explode? for k = 1:4 local f @@ -210,7 +240,8 @@ end @testset "Special solver solver_lsmr_se with Standard Error" begin se_solver = solver = (x, y) -> Unfold.solver_default(x, y, stderror = true) - m_tul_se = coeftable(fit(UnfoldModel, f, evts, data_r, basisfunction, solver = se_solver)) #with + m_tul_se = + coeftable(fit(UnfoldModel, f, evts, data_r, basisfunction, solver = se_solver)) #with m_tul = coeftable(fit(UnfoldModel, f, evts, data_r, basisfunction)) #without @test all(m_tul_se.estimate .== m_tul.estimate) @test !all(isnothing.(m_tul_se.stderror)) @@ -268,7 +299,7 @@ end if 1 == 0 # Fit mass-univariate 1st level for all subjects - basisfunction = firbasis(τ = (-.1, 0.4), sfreq = 10, name = "A") + basisfunction = firbasis(τ = (-0.1, 0.4), sfreq = 10, name = "A") resAll = DataFrame() f = @formula 0 ~ 1 + condA + condB for s in unique(evts.subject) diff --git a/test/fit_LMM.jl b/test/fit_LMM.jl index bb5cf4e2..65b9f909 100644 --- a/test/fit_LMM.jl +++ b/test/fit_LMM.jl @@ -23,7 +23,8 @@ # cut the data into epochs # TODO This ignores subject bounds data_e, times = Unfold.epoch(data = data, tbl = evts, τ = (-1.0, 1.9), sfreq = 10) - data_missing_e, times = Unfold.epoch(data = data_missing, tbl = evts, τ = (-1.0, 1.9), sfreq = 10) + data_missing_e, times = + Unfold.epoch(data = data_missing, tbl = evts, τ = (-1.0, 1.9), sfreq = 10) evts_e, data_e = Unfold.dropMissingEpochs(copy(evts), data_e) evts_missing_e, data_missing_e = Unfold.dropMissingEpochs(copy(evts), data_missing_e) @@ -44,7 +45,7 @@ rtol = 0.1, ) - + # with missing @time m_mum = fit( UnfoldModel, @@ -93,7 +94,7 @@ ) - evts.subjectB = evts.subject; + evts.subjectB = evts.subject evts1 = evts[evts.condA.==0, :] evts2 = evts[evts.condA.==1, :] @@ -111,7 +112,7 @@ X1_lmm = designmatrix(UnfoldLinearMixedModel, f1_lmm, evts1, b1) X2_lmm = designmatrix(UnfoldLinearMixedModel, f2_lmm, evts2, b2) - r = fit(UnfoldLinearMixedModelContinuousTime, X1_lmm + X2_lmm, data); + r = fit(UnfoldLinearMixedModelContinuousTime, X1_lmm + X2_lmm, data) df = coeftable(r) @test isapprox( @@ -134,51 +135,98 @@ end ## Condense check for multi channel, multi @testset "LMM multi channel, multi basisfunction" begin - data,evts = loadtestdata("testCase3", dataPath = (@__DIR__) * "/data") - transform!(evts,:subject=>categorical=>:subject); - data = vcat(data',data') - - bA0 = firbasis(τ=(-0.0,0.1),sfreq=10,name="bA0") - bA1 = firbasis(τ=(0.1,0.2),sfreq=10,name="bA1") - evts.subject2 = evts.subject - fA0 = @formula 0~1+condB + zerocorr(1|subject) - fA1 =@formula 0~1+condB + zerocorr(1|subject2) - m = fit(UnfoldModel, - Dict(0=>(fA0,bA0), - 1=>(fA1,bA1)), - evts,data,eventcolumn="condA") - - res = coeftable(m) - - @test all(last(.!isnothing.(res.group),8)) - @test all(last(res.coefname,8).=="(Intercept)") + data, evts = loadtestdata("testCase3", dataPath = (@__DIR__) * "/data") + transform!(evts, :subject => categorical => :subject) + data = vcat(data', data') + + bA0 = firbasis(τ = (-0.0, 0.1), sfreq = 10, name = "bA0") + bA1 = firbasis(τ = (0.1, 0.2), sfreq = 10, name = "bA1") + evts.subject2 = evts.subject + fA0 = @formula 0 ~ 1 + condB + zerocorr(1 | subject) + fA1 = @formula 0 ~ 1 + condB + zerocorr(1 | subject2) + m = fit( + UnfoldModel, + Dict(0 => (fA0, bA0), 1 => (fA1, bA1)), + evts, + data, + eventcolumn = "condA", + ) + + res = coeftable(m) + + @test all(last(.!isnothing.(res.group), 8)) + @test all(last(res.coefname, 8) .== "(Intercept)") end @testset "LMM bug reorder #115" begin - - data,evts = UnfoldSim.predef_2x2(;return_epoched=true,n_subjects=10,noiselevel=1) - -designList = [Dict(Any=>(@formula(0~1+A+B+zerocorr(1+B+A|subject)+zerocorr(1+B|item)),range(0,1,length=size(data,1)))), - Dict(Any=>(@formula(0~1+A+B+zerocorr(1+A+B|subject)+zerocorr(1+B|item)),range(0,1,length=size(data,1)))), - Dict(Any=>(@formula(0~1+zerocorr(1+A+B|subject)+zerocorr(1|item)),range(0,1,length=size(data,1))))] -#des = designList[1] -#des = designList[2] - for des = designList - @test_throws AssertionError fit(UnfoldModel,des,evts,data) + + data, evts = + UnfoldSim.predef_2x2(; return_epoched = true, n_subjects = 10, noiselevel = 1) + + designList = [ + Dict( + Any => ( + @formula( + 0 ~ + 1 + A + B + zerocorr(1 + B + A | subject) + zerocorr(1 + B | item) + ), + range(0, 1, length = size(data, 1)), + ), + ), + Dict( + Any => ( + @formula( + 0 ~ + 1 + A + B + zerocorr(1 + A + B | subject) + zerocorr(1 + B | item) + ), + range(0, 1, length = size(data, 1)), + ), + ), + Dict( + Any => ( + @formula(0 ~ 1 + zerocorr(1 + A + B | subject) + zerocorr(1 | item)), + range(0, 1, length = size(data, 1)), + ), + ), + ] + #des = designList[1] + #des = designList[2] + for des in designList + @test_throws AssertionError fit(UnfoldModel, des, evts, data) # end #counter check - des = Dict(Any=>(@formula(0~1+zerocorr(1|item)+zerocorr(1+A+B|subject)),range(0,1,length=size(data,1)))) - uf = fit(UnfoldModel,des,evts,data) - @test 3 == unique(@subset(coeftable(uf),@byrow(:group == Symbol("subject")),@byrow :time == 0.0).coefname) |> length + des = Dict( + Any => ( + @formula(0 ~ 1 + zerocorr(1 | item) + zerocorr(1 + A + B | subject)), + range(0, 1, length = size(data, 1)), + ), + ) + uf = fit(UnfoldModel, des, evts, data) + @test 3 == + unique( + @subset( + coeftable(uf), + @byrow(:group == Symbol("subject")), + @byrow :time == 0.0 + ).coefname, + ) |> length end @testset "LMM bug reshape #110" begin - data,evts = UnfoldSim.predef_2x2(;return_epoched=true,n_subjects=10,noiselevel=1) - des = Dict(Any=>(@formula(0~1+A+B+zerocorr(1+B+A|item)+zerocorr(1+B|subject)),range(0,1,length=size(data,1)))) - uf = fit(UnfoldModel,des,evts,data) - @test size(coef(uf)) ==(1,100,3) -end \ No newline at end of file + data, evts = + UnfoldSim.predef_2x2(; return_epoched = true, n_subjects = 10, noiselevel = 1) + des = Dict( + Any => ( + @formula( + 0 ~ 1 + A + B + zerocorr(1 + B + A | item) + zerocorr(1 + B | subject) + ), + range(0, 1, length = size(data, 1)), + ), + ) + uf = fit(UnfoldModel, des, evts, data) + @test size(coef(uf)) == (1, 100, 3) +end diff --git a/test/io.jl b/test/io.jl index 62f16599..f409e5c3 100644 --- a/test/io.jl +++ b/test/io.jl @@ -1,6 +1,6 @@ using UnfoldSim using DataFrames -save_path = mktempdir(;cleanup=false)#tempdir() +save_path = mktempdir(; cleanup = false)#tempdir() ## 1. Test data set: # - Generate a P1/N1/P3 complex for one subject (using UnfoldSim) @@ -9,84 +9,96 @@ save_path = mktempdir(;cleanup=false)#tempdir() -data1, evts1 = UnfoldSim.predef_eeg(; n_repeats=100, noiselevel=0.8); +data1, evts1 = UnfoldSim.predef_eeg(; n_repeats = 100, noiselevel = 0.8); # Assume that the data is from two different events evts1[!, :type] = repeat(["event_A", "event_B"], nrow(evts1) ÷ 2); -bf1_A = firbasis(τ=[-0.1, 1], sfreq=100, name="event_A"); -bf1_B = firbasis(τ=[-0.1, 1], sfreq=100, name="event_B"); +bf1_A = firbasis(τ = [-0.1, 1], sfreq = 100, name = "event_A"); +bf1_B = firbasis(τ = [-0.1, 1], sfreq = 100, name = "event_B"); f1_A = @formula 0 ~ 1; f1_B = @formula 0 ~ 1 + condition + spl(continuous, 4); bfDict1 = Dict("event_A" => (f1_A, bf1_A), "event_B" => (f1_B, bf1_B)); -data1_e,times = Unfold.epoch(data1,evts1,[-0.1,1],100) +data1_e, times = Unfold.epoch(data1, evts1, [-0.1, 1], 100) bfDict1_e = Dict("event_A" => (f1_A, times), "event_B" => (f1_B, times)); -for deconv = [false,true] +for deconv in [false, true] if deconv - - m1 = Unfold.fit(UnfoldModel, bfDict1, evts1, data1, eventcolumn="type"); + + m1 = Unfold.fit(UnfoldModel, bfDict1, evts1, data1, eventcolumn = "type") else - m1 = Unfold.fit(UnfoldLinearModel, bfDict1_e, evts1, data1_e; eventcolumn="type"); + m1 = Unfold.fit(UnfoldLinearModel, bfDict1_e, evts1, data1_e; eventcolumn = "type") end @testset "SingleSubjectDesign with two event types and splines" begin - # save the model to a compressed .jld2 file and load it again - - save(joinpath(save_path, "m1_compressed.jld2"), m1; compress=true) - m1_loaded = load(joinpath(save_path, "m1_compressed.jld2"), UnfoldModel, generate_Xs=true) - - @test ismissing(m1_loaded.designmatrix.Xs) == false - @test typeof(m1) == typeof(m1_loaded) - @test coeftable(m1) == coeftable(m1_loaded) - @test m1.modelfit.estimate == m1_loaded.modelfit.estimate - @test m1.designmatrix.events == m1_loaded.designmatrix.events - - # In the loaded version one gets two matrices instead of one. - # The dimensions do also not match. - # Probably the designmatrix reconstruction in the load function needs to be changed. - @test m1.designmatrix.Xs == m1_loaded.designmatrix.Xs - - # Test whether the effects function works with the loaded model - # and the results match the ones of the original model - eff1 = effects(Dict(:condition => ["car", "face"], :continuous => -5:1), m1) - eff1_loaded = effects(Dict(:condition => ["car", "face"], :continuous => -5:1), m1_loaded) - @test eff1 == eff1_loaded + # save the model to a compressed .jld2 file and load it again - # load the model without reconstructing the designmatrix - m1_loaded_without_dm = load(joinpath(save_path, "m1_compressed.jld2"), UnfoldModel, generate_Xs=false) + save(joinpath(save_path, "m1_compressed.jld2"), m1; compress = true) + m1_loaded = + load(joinpath(save_path, "m1_compressed.jld2"), UnfoldModel, generate_Xs = true) - - # ismissing should only be true fr the deconv case - @test ismissing(m1_loaded_without_dm.designmatrix.Xs) == (deconv == true) - -end + @test ismissing(m1_loaded.designmatrix.Xs) == false + @test typeof(m1) == typeof(m1_loaded) + @test coeftable(m1) == coeftable(m1_loaded) + @test m1.modelfit.estimate == m1_loaded.modelfit.estimate + @test m1.designmatrix.events == m1_loaded.designmatrix.events + + # In the loaded version one gets two matrices instead of one. + # The dimensions do also not match. + # Probably the designmatrix reconstruction in the load function needs to be changed. + @test m1.designmatrix.Xs == m1_loaded.designmatrix.Xs + + # Test whether the effects function works with the loaded model + # and the results match the ones of the original model + eff1 = effects(Dict(:condition => ["car", "face"], :continuous => -5:1), m1) + eff1_loaded = + effects(Dict(:condition => ["car", "face"], :continuous => -5:1), m1_loaded) + @test eff1 == eff1_loaded + + # load the model without reconstructing the designmatrix + m1_loaded_without_dm = load( + joinpath(save_path, "m1_compressed.jld2"), + UnfoldModel, + generate_Xs = false, + ) + + + # ismissing should only be true fr the deconv case + @test ismissing(m1_loaded_without_dm.designmatrix.Xs) == (deconv == true) + + end end #---- ## 2. Test data set: # - Generate a 2x2 design with Hanning window for multiple subjects (using UnfoldSim) # - Use a Mixed-effects Unfold model -data2, evts2 = UnfoldSim.predef_2x2(; n_subjects=5, return_epoched=true); +data2, evts2 = UnfoldSim.predef_2x2(; n_subjects = 5, return_epoched = true); # Define a model formula with interaction term and random effects for subjects f2 = @formula(0 ~ 1 + A * B + (1 | subject)); τ2 = [-0.1, 1]; sfreq2 = 100; -times2 = range(τ2[1], length=size(data2, 1), step=1 ./ sfreq2); +times2 = range(τ2[1], length = size(data2, 1), step = 1 ./ sfreq2); -m2 = Unfold.fit(UnfoldModel, Dict(Any => (f2, times2)), evts2, reshape(data2, 1, size(data2)...)); +m2 = Unfold.fit( + UnfoldModel, + Dict(Any => (f2, times2)), + evts2, + reshape(data2, 1, size(data2)...), +); -save(joinpath(save_path, "m2_compressed.jld2"), m2; compress=true) -m2_loaded = load(joinpath(save_path, "m2_compressed.jld2"), UnfoldModel, generate_Xs=true) +save(joinpath(save_path, "m2_compressed.jld2"), m2; compress = true) +m2_loaded = load(joinpath(save_path, "m2_compressed.jld2"), UnfoldModel, generate_Xs = true) @testset "2x2 MultiSubjectDesign Mixed-effects model" begin # save the model to a compressed .jld2 file and load it again - save(joinpath(save_path, "m2_compressed.jld2"), m2; compress=true) - m2_loaded = load(joinpath(save_path, "m2_compressed.jld2"), UnfoldModel, generate_Xs=true) + save(joinpath(save_path, "m2_compressed.jld2"), m2; compress = true) + m2_loaded = + load(joinpath(save_path, "m2_compressed.jld2"), UnfoldModel, generate_Xs = true) @test ismissing(m2_loaded.designmatrix.Xs) == false @test typeof(m2) == typeof(m2_loaded) @@ -106,7 +118,8 @@ m2_loaded = load(joinpath(save_path, "m2_compressed.jld2"), UnfoldModel, generat @test_broken eff2 == eff2_loaded # load the model without reconstructing the designmatrix - m2_loaded_without_dm = load(joinpath(save_path, "m2_compressed.jld2"), UnfoldModel, generate_Xs=false) + m2_loaded_without_dm = + load(joinpath(save_path, "m2_compressed.jld2"), UnfoldModel, generate_Xs = false) @test ismissing(m2_loaded_without_dm.designmatrix.Xs) == true -end \ No newline at end of file +end diff --git a/test/predict.jl b/test/predict.jl index e28ae413..f99b8f06 100644 --- a/test/predict.jl +++ b/test/predict.jl @@ -69,13 +69,19 @@ data, evts = loadtestdata("test_case_4a") # b1 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisA") b2 = firbasis(τ = (0.0, 0.95), sfreq = 20, name = "basisB") f = @formula 0 ~ 1 # 1 -m_tul = fit(UnfoldModel, Dict("eventA"=>(f,b1),"eventB"=>(f,b2)), evts, data,eventcolumn="type") +m_tul = fit( + UnfoldModel, + Dict("eventA" => (f, b1), "eventB" => (f, b2)), + evts, + data, + eventcolumn = "type", +) -p = predict(m_tul,DataFrame(:Cond => [1])) +p = predict(m_tul, DataFrame(:Cond => [1])) -@test size(p,1) == 40 -@test length(unique(p.time)) ==20 -@test unique(p.basisname) == ["basisA","basisB"] +@test size(p, 1) == 40 +@test length(unique(p.time)) == 20 +@test unique(p.basisname) == ["basisA", "basisB"] diff --git a/test/runtests.jl b/test/runtests.jl index 13b3fc01..51162dbc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,7 +28,7 @@ end include("splines.jl") end -@testset "Predict" begin +@testset "Predict" begin include("predict.jl") end diff --git a/test/setup.jl b/test/setup.jl index 233f43c2..cb98dfcf 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -8,8 +8,8 @@ using Random using CategoricalArrays using StatsBase -using MixedModels, RobustModels,BSplineKit # extensionTriggers +using MixedModels, RobustModels, BSplineKit # extensionTriggers using SparseArrays using UnfoldSim -include("test_utilities.jl") \ No newline at end of file +include("test_utilities.jl") diff --git a/test/splines.jl b/test/splines.jl index 55191ce7..99ef4af1 100644 --- a/test/splines.jl +++ b/test/splines.jl @@ -4,7 +4,7 @@ data, evts = loadtestdata("test_case_3a") # f_spl = @formula 0 ~ 1 + conditionA + spl(continuousA, 4) # 1 f = @formula 0 ~ 1 + conditionA + continuousA # 1 data_r = reshape(data, (1, :)) -data_e, times = Unfold.epoch(data=data_r, tbl=evts, τ=(-1.0, 1.0), sfreq=10) # cut the data into epochs +data_e, times = Unfold.epoch(data = data_r, tbl = evts, τ = (-1.0, 1.0), sfreq = 10) # cut the data into epochs m_mul = coeftable(fit(UnfoldModel, f, evts, data_e, times)) m_mul_spl = coeftable(fit(UnfoldModel, f_spl, evts, data_e, times)) @@ -20,12 +20,12 @@ s = Unfold.formula(fit(UnfoldModel, f_spl, evts, data_e, times)).rhs.terms[3] @testset "outside bounds" begin # test safe prediction m = fit(UnfoldModel, f_spl, evts, data_e, times) - r = predict(m, DataFrame(conditionA=[0, 0], continuousA=[0.9, 1.9])) + r = predict(m, DataFrame(conditionA = [0, 0], continuousA = [0.9, 1.9])) @test all(ismissing.(r.yhat[r.continuousA.==1.9])) @test !any(ismissing.(r.yhat[r.continuousA.==0.9])) end -basisfunction = firbasis(τ=(-1, 1), sfreq=10, name="A") +basisfunction = firbasis(τ = (-1, 1), sfreq = 10, name = "A") @testset "timeexpanded" begin # test time expanded m_tul = coeftable(fit(UnfoldModel, f, evts, data_r, basisfunction)) @@ -35,7 +35,7 @@ end # test safe predict m = fit(UnfoldModel, f_spl, evts, data_r, basisfunction) - p = predict(m, DataFrame(conditionA=[0, 0, 0], continuousA=[0.9, 0.9, 1.9])) + p = predict(m, DataFrame(conditionA = [0, 0, 0], continuousA = [0.9, 0.9, 1.9])) @test all(ismissing.(p[p.continuousA.==1.9, :yhat])) end @@ -49,7 +49,7 @@ if 1 == 0 using AlgebraOfGraphics yhat_mul.conditionA = categorical(yhat_mul.conditionA) yhat_mul.continuousA = categorical(yhat_mul.continuousA) - m = mapping(:times, :yhat, color=:continuousA, linestyle=:conditionA) + m = mapping(:times, :yhat, color = :continuousA, linestyle = :conditionA) df = yhat_mul AlgebraOfGraphics.data(df) * visual(Lines) * m |> draw end @@ -66,7 +66,7 @@ end f_evaluated = Unfold.formula(m) effValues = [-1, -0.99, 0, 0.99, 1] - effValues = range(-1.1, 1.1, step=0.1) + effValues = range(-1.1, 1.1, step = 0.1) effSingle = effects(Dict(:continuousA => effValues), m) tmp = subset(effSingle, :time => x -> x .== -1.0) @test tmp.yhat[tmp.continuousA.==-1.1] ≈ tmp.yhat[tmp.continuousA.==0.9] @@ -79,4 +79,4 @@ end f_spl = @formula 0 ~ 1 + conditionA + spl(continuousA, 3) # 1 @test_throws AssertionError fit(UnfoldModel, f_spl, evts, data_e, times) -end \ No newline at end of file +end diff --git a/test/statistics.jl b/test/statistics.jl index ddd04db8..827d05b4 100644 --- a/test/statistics.jl +++ b/test/statistics.jl @@ -1,33 +1,38 @@ @testset "LMM LRT" begin - data,evts = loadtestdata("testCase3",dataPath="data"); - data = reshape(data,1,length(data)) # add second channel - data = vcat(data,data); - - evts[!,:subject] .= string.(evts.subject); - data_epoch,times = Unfold.epoch(data=data.+1.1.*randn(MersenneTwister(1),size(data)),tbl=evts,τ=(-0.1,0.3),sfreq=10); - evts_epoch,data_epoch = Unfold.dropMissingEpochs(evts,data_epoch) - - - f0 = @formula 0~1+condA + (1|subject); - f1 = @formula 0~1+condA+condB + (1|subject); - - m0 = fit(UnfoldModel,Dict(Any=>(f0,times)),evts_epoch,data_epoch); - m1 = fit(UnfoldModel,Dict(Any=>(f1,times)),evts_epoch,data_epoch); - - evts[!,:y] = data_epoch[1,1,:] - - f0 = @formula y~1+condA + (1|subject); - f1 = @formula y~1+condA+condB + (1|subject); - - lmm0 = fit(MixedModel,f0,evts) - lmm1 = fit(MixedModel,f1,evts) - - uf_lrt = likelihoodratiotest(m0,m1) - mm_lrt = MixedModels.likelihoodratiotest(lmm0,lmm1) - - @test mm_lrt.pvalues ≈ uf_lrt[1].pvalues - - - - -end \ No newline at end of file + data, evts = loadtestdata("testCase3", dataPath = "data") + data = reshape(data, 1, length(data)) # add second channel + data = vcat(data, data) + + evts[!, :subject] .= string.(evts.subject) + data_epoch, times = Unfold.epoch( + data = data .+ 1.1 .* randn(MersenneTwister(1), size(data)), + tbl = evts, + τ = (-0.1, 0.3), + sfreq = 10, + ) + evts_epoch, data_epoch = Unfold.dropMissingEpochs(evts, data_epoch) + + + f0 = @formula 0 ~ 1 + condA + (1 | subject) + f1 = @formula 0 ~ 1 + condA + condB + (1 | subject) + + m0 = fit(UnfoldModel, Dict(Any => (f0, times)), evts_epoch, data_epoch) + m1 = fit(UnfoldModel, Dict(Any => (f1, times)), evts_epoch, data_epoch) + + evts[!, :y] = data_epoch[1, 1, :] + + f0 = @formula y ~ 1 + condA + (1 | subject) + f1 = @formula y ~ 1 + condA + condB + (1 | subject) + + lmm0 = fit(MixedModel, f0, evts) + lmm1 = fit(MixedModel, f1, evts) + + uf_lrt = likelihoodratiotest(m0, m1) + mm_lrt = MixedModels.likelihoodratiotest(lmm0, lmm1) + + @test mm_lrt.pvalues ≈ uf_lrt[1].pvalues + + + + +end diff --git a/test/time_expand.jl b/test/time_expand.jl index df9fda8d..51adf394 100644 --- a/test/time_expand.jl +++ b/test/time_expand.jl @@ -8,7 +8,7 @@ term = Unfold.TimeExpandedTerm(FormulaTerm(Term, Term), basisfunction, :latency) Xdc = Unfold.time_expand(X, term, tbl) kernel = Unfold.kernel -ncolsBasis = size(kernel(term.basisfunction)(0.), 2) +ncolsBasis = size(kernel(term.basisfunction)(0.0), 2) X = reshape(X, size(X, 1), :) ncolsX = size(X)[2] @@ -22,8 +22,9 @@ else bases = kernel(term.basisfunction).(eachrow(tbl[!, term.eventfields])) end -rows = Unfold.timeexpand_rows(onsets,bases,Unfold.shiftOnset(term.basisfunction),ncolsX) -cols = Unfold.timeexpand_cols(term,bases,ncolsBasis,ncolsX) -vals = Unfold.timeexpand_vals(bases,X,size(cols),ncolsX) +rows = Unfold.timeexpand_rows(onsets, bases, Unfold.shiftOnset(term.basisfunction), ncolsX) +cols = Unfold.timeexpand_cols(term, bases, ncolsBasis, ncolsX) +vals = Unfold.timeexpand_vals(bases, X, size(cols), ncolsX) -@test Unfold.timeexpand_cols_allsamecols(bases,ncolsBasis,ncolsX) == Unfold.timeexpand_cols_generic(bases,ncolsBasis,ncolsX) +@test Unfold.timeexpand_cols_allsamecols(bases, ncolsBasis, ncolsX) == + Unfold.timeexpand_cols_generic(bases, ncolsBasis, ncolsX) diff --git a/test/utilities.jl b/test/utilities.jl index a42041c0..a56a6d74 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -1,33 +1,31 @@ @testset "epoch" begin d = collect(1:100) - evt = DataFrame(:latency=>(50)) - - ep = τ -> Unfold.epoch(d,evt,τ,1)[1][1,:,1] + evt = DataFrame(:latency => (50)) + + ep = τ -> Unfold.epoch(d, evt, τ, 1)[1][1, :, 1] # check time delays - @test ep((0,10.)) ≈ collect(50:60.) - @test ep((-10,10.)) ≈ collect(40:60.) - @test ep((-10,0.)) ≈ collect(40:50.) - @test ep((5,15)) ≈ collect(55:65.) - @test ep((-15,-5)) ≈ collect(35:45.) + @test ep((0, 10.0)) ≈ collect(50:60.0) + @test ep((-10, 10.0)) ≈ collect(40:60.0) + @test ep((-10, 0.0)) ≈ collect(40:50.0) + @test ep((5, 15)) ≈ collect(55:65.0) + @test ep((-15, -5)) ≈ collect(35:45.0) # check corner cases (sample doesnt end on sampling rate) - @test ep((0.6,2)) ≈ collect(51:52.) - @test ep((0.2,2)) ≈ collect(50:52.) + @test ep((0.6, 2)) ≈ collect(51:52.0) + @test ep((0.2, 2)) ≈ collect(50:52.0) # test sampling frequencies - ep = τ -> Unfold.epoch(d,evt,τ,2)[1][1,:,1] - @test ep((-1.0, 2)) ≈ collect(48:54.) + ep = τ -> Unfold.epoch(d, evt, τ, 2)[1][1, :, 1] + @test ep((-1.0, 2)) ≈ collect(48:54.0) - ep = τ -> Unfold.epoch(d,evt,τ,0.5)[1][1,:,1] - @test ep((-4.0, 8)) ≈ collect(48:54.) + ep = τ -> Unfold.epoch(d, evt, τ, 0.5)[1][1, :, 1] + @test ep((-4.0, 8)) ≈ collect(48:54.0) -# rounding bug when latency was .5 -> bug #78 -d = zeros((1, 1270528)) -evt = DataFrame(:latency=>(181603.5)) -ep = τ -> Unfold.epoch(d,evt,τ,256.)[1][1,:,1] -ep((-0.1, 0.8)) + # rounding bug when latency was .5 -> bug #78 + d = zeros((1, 1270528)) + evt = DataFrame(:latency => (181603.5)) + ep = τ -> Unfold.epoch(d, evt, τ, 256.0)[1][1, :, 1] + ep((-0.1, 0.8)) end - -