The
file Family.txt
is from chapter 2 of Demidenko's book. Among other things, it
gives the heights and weights of 71 persons. A simple plot will
show a scattered increasing trend.
dat = read.table("Family.txt",head=T)
names(dat)
attach(dat)
plot(Weight,Height)
We can, of course, fit a regression to explore the height-weight
relationship:
$$
w_i = \alpha + \beta h_i + \epsilon_i.
$$
But here we have another relevant piece of info:
each of the 71 persons belongs to one of 18 families, and the
data set gives us the FamilyID's. Once would expect that the
members of the same family would behave similarly. Can we somehow
incorporate this information by adding a "family effect" to the model?
The first approach could be to fit a model:
$$
w_{ki} = \alpha_k + \beta_k h_{ki} + \epsilon_{ki}.
$$
Here each person is coded as a pair $(k,i)$,
meaning $i$-th person from the $k$-th family.
But this is not good because of the following reasons.
It is difficult to use the estimated parameters. For
instance, $\h \alpha_1$ is specific for family 1. And why
should one
care about that particular family?
Most families have a very small number of members in the
data. So the estimates of $\alpha_k$'s and $\beta_k$'s
will be bad (or impossible, in case of family 16 that has ony
one representative!).
A better approach is to allow the intercept to be family-specific
and random:
$$
\alpha_k = \mu + a_k,
$$
where $\mu$ is fixed and unknown and $a_k$'s are
iid $N(0,\sigma^2_a).$ So the full model now becomes
$$
w_{ki} = \mu+a_k + \beta h_{ki} + \epsilon_{ki}.
$$
If we club the random terms in the RHS together, then this is
actually a GLS model where the covariance matrix is a block
diagonal matrix with the one block per family. The $i$-th
block has $\sigma^2+\sigma^2_a$ in the diagonal entries and
the only $\sigma^2_a$ in the off-diagonal ones.
We can fit this model using R as follows.
library(nlme) #install it first if you have not done so already
fit=lme (Weight~Height , random=~ 1 |FamilyID)
summary(fit)
This example is based on the online Netmaster Statistics Courses. Here we
want to compare
two methods (HPLC and NIR) to ascertain the amount of active content
in tablets. Suppose that we want to test if the two methods yield more or less the
same measurement.
The tests have been applied to the same set of 10 tablets (e.g.,
breaking each tablet into two halves, and applying one method to each
half). The resulting data are shown in the following table.
Tablet $(i)$
HPLC $(x_i)$
NIR$(y_i)$
1
10.4
10.1
2
10.6
10.8
3
10.2
10.2
4
10.1
9.9
5
10.3
11.0
6
10.7
10.5
7
10.3
10.2
8
10.9
10.9
9
10.1
10.4
10
9.8
9.9
One standard method to analyze the data is to perform a paired
$t$-test, which basically works with the differences $z_i = x_i-y_i.$
dat = read.table('tablet.txt',head=T)
with(dat,t.test(HPLC,NIR,paired=T))
Now here it is quite natural to assume that $x_i$'s
and $y_i$'s will be positively correlated. This information
is not used in the paired $t$-test. However, we may use a
linear mixed model to incorporate the info. We shall allow a
tablet effect.
$$
y_{ij} = \mu + a_i + \beta_j + \epsilon_{ij}.
$$
Here $y_{ij}$ is the measurement for tablet $i$ under
method $j.$ The tablet effect $a_i$ is random. To fit
this model using R we need to rearrange the data:
How can you use LME to estimate a fixed
unknown length which has been measured 5 times each by 3
persons?
We are inspecting the quality of a survey instrument for measuring angles. 50
roughly equilaterial triangles are taken and each angle is measured
independently. We are testing
$H_0:$ sum of all the angles
of the same triangle equals $180^\circ.$
A straight-forward method is to apply $t$-test to the total
measured angles for the triangles. Can you improve upon this by
using LME?
Comments
To post an anonymous comment, click on the "Name" field. This
will bring up an option saying "I'd rather post as a guest."