# Accounting for groups in GzLM Poisson regression

I am currently trying to perform a regression on data with one predictor and one covariate.

The hypothesis I am trying to test is whether the incidence of pores (a hole through a cell) increases when I stretch the cells. I have three cell lines and three stretch levels (0, 10 and 20%). At each stretch level and for each cell line I have 4 (and sometimes 3).

Because I am sampling the occurrence of rare events (typical counts between 0 and 20) over an area, we are looking at a Poisson process. This is the reason why I assume the data is Poisson-distributed and not normally distributed. Therefore I am turning towards Poisson regression, a kind of Generalized linear modeling (GzLM).

x1=1x34 vector containing the strain-level for each specimen x2=1x34 vector containing the cell line fore each specimen y=1x34 vector containing the pore counts for each specimen There are 3x3x4-2=34 specimens: cell line 1 is missing one specimen at 10%, cell line 2 is missing one specimen at 20%.

I am doing my statistical analysis in MATLAB. Before I came to the conclusion that my data is Poisson distributed, I used ANCOVA with the aoctool function, fitting separate lines through each cell line. The aoctool then calculates an F-statistic and a corresponding p-value.

[h,atab,ctab,astats]=aoctool2(x1,yT,x2,0.05,'Strain','Pore Count','Cell Line','off','separate lines');


The aoctool fits the data to the following model: y = (α + αi) + (β + βi)x + ε where y is the dependent variable, α the general intercept, αi the cell line specific intercept, β the general slope, βi the cell line specific intercept and ε the error term.

Fitting the data to this model regresses the pore count with the stretch-level whilst accounting for possible differences between cell lines (groups).

Now I want to do the same with the fitglm (or glmfit) function in Matlab but using a Poisson rather than a normal distribution. However, I have not been able to get the fitglm function to fit separate lines for each cell line. Rather, it takes the cell line as a covariate.

My question is twofold:

1. Is there a difference between using cell line as a covariate or as grouping variable?
2. How do I adjust the fitglm-input to take into account cell line as a grouping variable?

This is what I have tried so far:

tbl=table(yT,x1,x2); % enter my data into a 'table' format for the fitglm function
tbl.x2=nominal(tbl.x2); % x2, cell line, is a nominal variable with values 1, 2 or 3
model=fitglm(tbl,'yT ~ (x1+x2)','Distribution','Poisson')


this results in

Estimated Coefficients: Estimate SE tStat pValue
________ _________ _______ __________

(Intercept)      2.5086     0.097237     25.798    9.234e-147
x1             0.024509    0.0024828     9.8714    5.5397e-23
x2_2           -0.47372      0.10328    -4.5869     4.499e-06
x2_3           -0.65656      0.10444    -6.2862    3.2532e-10


I seem to have found a solution. The grouping variable must be turned into a categorical variable and the model 'Pore Count ~ Group+Strain+Group*Strain' is equivalent to the separate lines fit from aoctool. Removing the interaction term (Group*Strain) will give a model equivalent to parallel lines fit in aoctool. The p-value on your interaction thus tells you whether the slopes of your separate lines are indeed significantly different.

Here is the basic model fit for the Poisson model:

ds=dataset(categorical(data.Eye),data.StdTI,data.(i+2),'VarNames',{'Group','Covariate','Response'});
gzlm=fitglm(ds,'Response ~ Group+Covariate+Group*Covariate','Distribution','poisson');


Here is code to show how fitlm (not fitglm) and aoctool give you the same result:

% I have been using the 'aoctool' function to perform ANCOVA with separate
% lines for each group and have been trying to calculate the R-squared value
% on the model (R2=1-SSresidual/SStotal). In doing some searches online I
% realised I could do the same fit with the function 'fitlm'. When I compare
% the results from aoctool and fitlm functions I get different results in the
% sum of squares of the model components (intercept, slope, interaction). The
% regressions are exactly the same, as is the residual sum of squares. The
% fitlm results match exactly with an SPSS-ANCOVA analysis I did. Does anyone
% know where this discrepancy comes from? Different type of sum of squares or
% is aoctool not adjusted/corrected, whereas fitlm is?

% ANSWER: use the component table of the ANOVA function and the SS are
% equivalent. Now it is time to calculate the R2

clear all
close all
COV=[1:6,2:7,3:8]'; %Covariate
Response=[COV(1:6)'.*normrnd(4,0.3,1,6),COV(7:12)'.*normrnd(2,0.2,1,6),...
COV(13:18)'.*normrnd(3,0.5,1,6)]'; %Response variable
Group=categorical([ones(1,6),2*ones(1,6),3*ones(1,6)]'); %Grouping
ds=dataset(Group,COV,Response,'VarNames',{'Group','Covariate','Response'});
lm=fitlm(ds,'Response ~ Group+Covariate+Group*Covariate');
[h,atab,ctab,stats]=aoctool(COV,Response,Group,0.05,...
'Covariate','Response','Group','on','separate lines');
lmtable=anova(lm,'component');
% Now compare 'atab' and lmtable