Fit a multivariate regression model to panel data, assuming different intercepts and common slopes.

Load the sample data.

load('flu')

The dataset array flu contains national CDC flu estimates, and nine separate regional estimates based on Google® query data.

Extract the response and predictor data.

Y = double(flu(:,2:end-1));
[n,d] = size(Y);
x = flu.WtdILI;

The responses in Y are the nine regional flu estimates. Observations exist for every week over a one-year period, so $$n$$ = 52. The dimension of the responses corresponds to the regions, so $$d$$ = 9. The predictors in x are the weekly national flu estimates.

Plot the flu data, grouped by region.

figure;
regions = flu.Properties.VarNames(2:end-1);
plot(x,Y,'x')
legend(regions,'Location','NorthWest')

Fit the multivariate regression model $$y_{ij}={\alpha }_j+\beta x_{ij}+{ϵ}_{ij}$$, where $$i=1,\dots ,n$$ and $$j=1,\dots ,d$$, with between-region concurrent correlation $$COV({ϵ}_{ij},{ϵ}_{ij})={\sigma }_{jj}$$.

There are $$K$$ = 10 regression coefficients to estimate: nine intercept terms and a common slope. The input argument X should be an $$n$$-element cell array of $$d$$ -by- $$K$$ design matrices.

X = cell(n,1);
for i = 1:n
	X{i} = [eye(d) repmat(x(i),d,1)];
end
[beta,Sigma] = mvregress(X,Y);

beta contains estimates of the $$K$$-dimensional coefficient vector $$({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_9,\beta {)}^{\prime }$$.

Sigma contains estimates of the $$d$$ -by- $$d$$ variance-covariance matrix $${({\sigma }_{ij})}_{d\times d}$$, $$i,j=1,\dots ,d$$ for the between-region concurrent correlations.

Plot the fitted regression model.

B = [beta(1:d)';repmat(beta(end),1,d)];
xx = linspace(.5,3.5)';
fits = [ones(size(xx)),xx]*B;

figure;
h = plot(x,Y,'x',xx,fits,'-');
for i = 1:d
	set(h(d+i),'color',get(h(i),'color'));
end
legend(regions,'Location','NorthWest');

The plot shows that each regression line has a different intercept but the same slope. Upon visual inspection, some regression lines appear to fit the data better than others.

Fit a multivariate regression model to panel data using least squares, assuming different intercepts and slopes.

Load the sample data.

load('flu');

The dataset array flu contains national CDC flu estimates, and nine separate regional estimates based on Google® queries.

Extract the response and predictor data.

Y = double(flu(:,2:end-1));
[n,d] = size(Y);
x = flu.WtdILI;

The responses in Y are the nine regional flu estimates. Observations exist for every week over a one-year period, so $$n$$ = 52. The dimension of the responses corresponds to the regions, so $$d$$ = 9. The predictors in x are the weekly national flu estimates.

Fit the multivariate regression model $$y_{ij}={\alpha }_j+{\beta }_j{x}_{ij}+{ϵ}_{ij}$$, where $$i=1,\dots ,n$$ and $$j=1,\dots ,d$$, with between-region concurrent correlation $$COV({ϵ}_{ij},{ϵ}_{ij})={\sigma }_{jj}$$.

There are $$K$$ = 18 regression coefficients to estimate: nine intercept terms, and nine slope terms. X is an $$n$$-element cell array of $$d$$ -by- $$K$$ design matrices.

X = cell(n,1);
for i = 1:n
    X{i} = [eye(d) x(i)*eye(d)];
end
[beta,Sigma] = mvregress(X,Y,'algorithm','cwls');

beta contains estimates of the $$K$$-dimensional coefficient vector $$({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_9,{\beta }_1,{\beta }_2,\dots ,{\beta }_9{)}^{\prime }$$.

Plot the fitted regression model.

B = [beta(1:d)';beta(d+1:end)'];
xx = linspace(.5,3.5)';
fits = [ones(size(xx)),xx]*B;

figure;
h = plot(x,Y,'x',xx,fits,'-');
for i = 1:d
	set(h(d+i),'color',get(h(i),'color'));
end

regions = flu.Properties.VarNames(2:end-1);
legend(regions,'Location','NorthWest');

The plot shows that each regression line has a different intercept and slope.

Fit a multivariate regression model using a single $$n$$-by-$$P$$ design matrix for all response dimensions.

Load the sample data.

load('flu')

The dataset array flu contains national CDC flu estimates, and nine separate regional estimates based on Google® queries.

Extract the response and predictor data.

Y = double(flu(:,2:end-1));
[n,d] = size(Y);
x = flu.WtdILI;

The responses in Y are the nine regional flu estimates. Observations exist for every week over a one-year period, so $$n$$ = 52. The dimension of the responses corresponds to the regions, so $$d$$ = 9. The predictors in x are the weekly national flu estimates.

Create an $$n$$ -by- $$P$$ design matrix X. Add a column of ones to include a constant term in the regression.

X = [ones(size(x)),x];

Fit the multivariate regression model

where $$i=1,\dots ,n$$ and $$j=1,\dots ,d$$, with between-region concurrent correlation

There are 18 regression coefficients to estimate: nine intercept terms, and nine slope terms.

[beta,Sigma,E,CovB,logL] = mvregress(X,Y);

beta contains estimates of the $$P$$-by-$$d$$ coefficient matrix. Sigma contains estimates of the $$d$$-by-$$d$$ variance-covariance matrix for the between-region concurrent correlations. E is a matrix of the residuals. CovB is the estimated variance-covariance matrix of the regression coefficients. logL is the value of the log likelihood objective function after the last iteration.

Plot the fitted regression model.

B = beta;
xx = linspace(.5,3.5)';
fits = [ones(size(xx)),xx]*B;

figure
h = plot(x,Y,'x', xx,fits,'-');
for i = 1:d
    set(h(d+i),'color',get(h(i),'color'))
end

regions = flu.Properties.VarNames(2:end-1);
legend(regions,'Location','NorthWest')

The plot shows that each regression line has a different intercept and slope.