Тренинг Платформа Engee (LTC/MBD)
Интеграция с другими языками программирования
Пример NLME: Реализация NLME модели в окружении Engee
В этом примере мы предлагаем вам запустить пример по оценке моделей с произвольными смешанными взаимосвязями (NLME, Nonlinear with Mixed Effects). Пример разработан для MATLAB и изложен в справочной системе по команде nlmefit.
using MATLAB
mat"ver -support"
Выполним пример nlmefit
Задать график роста нескольких саженцев апельсиновых деревьев.
out = mat"CIRC = [30 58 87 115 120 142 145;
33 69 111 156 172 203 203;
30 51 75 108 115 139 140;
32 62 112 167 179 209 214;
30 49 81 125 142 174 177];
time = [118 484 664 1004 1231 1372 1582];
[time;CIRC]
"
Как обычно, мы получаем на вывод в Engee именно ту строку, которая бы отобразилась после выполнения последней команды ячейки Live скрипта MATLAB.
1. Выведем график роста саженцев
В качестве примера можно построить график в MATLAB и сохранить график во внешний файл. В этом примере мы построим график средствами Julia:
using Plots
mtime = out[1,:];
CIRC = out[2:end,:]';
scatter(mtime, CIRC, leg=:topleft, markersize=8, markerstrokecolor=:white)
Зададим общий вид модели, которая будет предсказывать рост саженцев
mat"model = @(PHI,t)(PHI(:,1))./(1+exp(-(t-PHI(:,2))./PHI(:,3)));"
2. Подберем коэффициенты модели
В MATLAB эту операцию можно выполнить при помощи функции nlmefit:
mat"
TIME = repmat(time,5,1);
NUMS = repmat((1:5)',size(time));
beta0 = [100 100 100];
[beta1,PSI1,stats1] = nlmefit(TIME(:),CIRC(:),NUMS(:), [],model,beta0)
disp('Done')
"
3. Проанализируем результат и упростим модель
Как видно из вывода, влияние второго по счету предиктора (PSI1(2,2)) пренебрежимо мало. Устраним его из модели, задав аргумент REParamsSelect:
mat"
[beta2,PSI2,stats2,b2] = nlmefit(TIME(:),CIRC(:), NUMS(:),[],model,beta0,'REParamsSelect',[1 3])
disp('Done')
"
Параметр правдоподобия logl не изменился, AIC и BIC (критерии Акаике и Байеса) уменьшились, так что решение устранить второй предиктор было правильным.
4. Построим график
Теперь нам нужно рассчитать графики роста по модели beta2, в которую входит оценка случайных воздействий, и посмотреть, будет ли график близок к реальным данным:
pp = mat"
PHI = repmat(beta2,1,5) + ... % Fixed effects
[b2(1,:);zeros(1,5);b2(2,:)]; % Random effects
tplot = 0:1:1600;
pplot = zeros(length(tplot), 5);
for i = 1:5
fitted_model=@(t)(PHI(1,i))./(1+exp(-(t-PHI(2,i))./PHI(3,i)));
pplot(:,i) = fitted_model(tplot);
end";
Получим матрицу pplot из окружения MATLAB в окружение Julia:
pp = @mget pplot;
И построим наглядный график в среде Engee: