In this article, we’ll cover the main functions to conduct a statistical analysis. I’m working on a sample data set that describes the salaries of employees depending on several features.

# I. Data exploration

The file I’m working on comes with a small description DES. The description of the data is the following.

wage      educ      exper     tenure    nonwhite  female    married   numdep
smsa      northcen  south     west      construc  ndurman   trcommpu  trade
services  profserv  profocc   clerocc   servocc   lwage     expersq   tenursq

Obs:   526

1. wage                     average hourly earnings
2. educ                     years of education
3. exper                    years potential experience
4. tenure                   years with current employer
5. nonwhite                 =1 if nonwhite
6. female                   =1 if female
7. married                  =1 if married
8. numdep                   number of dependents
9. smsa                     =1 if live in SMSA
10. northcen                 =1 if live in north central U.S
11. south                    =1 if live in southern region
12. west                     =1 if live in western region
13. construc                 =1 if work in construc. indus.
14. ndurman                  =1 if in nondur. manuf. indus.
15. trcommpu                 =1 if in trans, commun, pub ut
16. trade                    =1 if in wholesale or retail
17. services                 =1 if in services indus.
18. profserv                 =1 if in prof. serv. indus.
19. profocc                  =1 if in profess. occupation
20. clerocc                  =1 if in clerical occupation
21. servocc                  =1 if in service occupation
22. lwage                    log(wage)
23. expersq                  exper^2
24. tenursq                  tenure^2


The first step is to load the data set in our path :

wage1=load('WAGE1.raw')
[n,k]=size(wage1)


Then, we’ll plot a histogram of hourly salaries :

hist(wage1(:,1))


Now, let’s plot the histogram for the number of years of education.

We can take a look at the mean of each feature :

mean(wage1(:,:))'

ans =

5.8961
12.5627
17.0171
5.1046
0.1027
...


We can also take a look at other descriptive statistics :

std(wage1(:,:))'
min(wage1)'
max(wage1)'
cov(wage1)'
corrcoef(wage1(:,:))'


# II. Linear regression

## Compute coefficients

We now want to assess the effect of columns 2, 3 and 4 (education, experience, years with current employer) on the salary.

$y_i = \beta_0 + \beta_1 * {x_1}_i + \beta_2 * {x_2}_i + \beta_3 * {x_3}_i + u_i$
y=wage1(:,1)
X=[ones(n,1),wage1(:,[2,3,4])]


We can compute the Beta coefficients and the residuals :

beta=inv(X'*X)*X'*y
u=y-X*beta


The coefficients are the following :

beta =

-2.8727
0.5990
0.0223
0.1693


The first coefficient is constant. Then, the model indicates that an additional year of education should increase the hourly salary by 0.60\$.

## Residuals

### Effect of each feature individually

We can also plot a histogram of the residuals. The distribution of the residual appears to be close to a normal one.

We can decide to remove outliers above 5 for example :

indices = find(u>2.5);
u(indices) = [];
hist(u)


## Test hypothesis

Are the coefficients significant?

To assess the significance of each estimator, we need to estimate their standard deviation.

sig2=u'*u/(n-4)
std=sqrt(diag(sig2*inv(X'*X)))


Then, we can easily compute the test statistic :

t = beta./std

t =

-5.6300
16.6857
2.6471
11.1725


What should we compare the test statistics? Well, to the critical value! For example, since we estimate 4 parameters, the 95% critical value is :

C1 = tinv(0.95, n-4)

C1 =

1.6478


The T-stats are all greater than the critical value. For this reason, we reject the hypothesis $$H_0$$ that the variables have no impact on the output.

### Effect of a single feature

If we are interested in the effect of a single feature, we can specify a value for its effect for example.

t_educ = (beta(2)-0.1)/std(2)


### Effect of each feature on log(wage)

Testing the effect of the features on the log(wage) is a way to assess their importance in terms of the percentage change.

y=log(wage1(:,1))
beta=inv(X'*X)*X'*y
std=sqrt(diag(sig2*inv(X'*X)))
t = beta./std


### Test is 2 features have the same effect

Now, let’s test the following hypothesis: the effect of education is the same as the effect of professional experience. This can be rewritten as : $$H_0 : \beta_2 = \beta_3$$ . However, estimating the variance of $$\beta_2 - \beta_3$$ to build the test hypothesis would be challenging. We can apply a simple trick by building a new feature as the sum of the 2 features to test.

First of all, append the new feature to the model :

capitaltot = X(:,2) + X(:,3)


This introduces some multicolinearity which we should control by removing $$X_3$$ :

X=[X(:,[1,2,4]),capitaltot];
y=log(wage1(:,1))


We can run the test : $$H_0 : \beta_4 = 0$$

beta=inv(X'*X)*X'*y
std=sqrt(diag(sig2*inv(X'*X)))
t = beta(4)/std(4)


t =

0.4883


The T-test is smaller than the critical value. For this reason, we cannot reject the hypothesis that the effect in the percentage of education is the same as the one of professional experience.

### Fisher Test for joint hypothesis

Let us now consider the hypothesis : $$H_0 : \beta_2 = 0, \beta_3 = 0$$ . In such case, we should not apply 2 independent tests, since we are interested in the overall hypothesis.

We need to define a non-constrained model, a constrained model and test the difference between those two models.

% non-constrained model
X0=[ones(n,1),wage1(:,[2,3,4])]
y0=log(wage1(:,1))
beta0=inv(X0'*X0)*X0'*y0
u0=y0-X0*beta0
SSR0=u0'*u0

% model contraint
X1=[ones(n,1),wage1(:,[4])]
y1=log(wage1(:,1))
beta1=inv(X1'*X1)*X1'*y1
u1=y1-X1*beta1
SSR1=u1'*u1


The Fisher F-Stat is defined as :

$\hat{F}_j = \frac { {SSR}_1 - {SSR}_0 } { {SSR}_0} \frac {n - K} {2}$

Where K is the number of variables included in the non-constrained model (4 here).

K = 4
F = (SSR1-SSR0)/(SSR0) * (n-K)/2


F =

80.1478


### Interaction effect

Let’s define new variables :

X=wage1
educ=X(:,2)
exper=X(:,3)
tenure=X(:,4)
female=X(:,6)
married=X(:,7)
marrmale = (1-female).*married
marrfem = female.*married
singfem=female.*(1-married)


We create interaction variables :

• effect of a married male
• effect of a married female
• effect of a single female

The base case is a single male. We can define test statistics on those variables as we did above.

**Conclusion **: This was a quick introduction to Matlab. I hope you found it useful. Don’t hesitate to drop a comment if you have a question.