[R-sig-ME] Simple one level model with known variances

[robert.h.creecy at census.gov]

All,

I have a very simple set of county data, with one dependent variable (y)
and three covariates (x),
one observation for each county. I have a lot of these type of problems to
do.

y is actually a summary statistic that comes from a survey within the
county,
and a variable (psi2) of the estimated sampling variance of y for each
county is also provided.

I am trying to figure out how to specify the known variances psi2_i within
the model

y_i = X_i beta + u_i + e_i    ,   b_i ~ N (0, psi2_i)  ,     e_i ~ N(0,
sigma^2)

u_i is a random effect for county i that comes from the uncertainty in the
survey estimate for that county.

I have looked at Chap. 5 of Pinheiro and Bates trying out how to use the
variance functions
as weights in lme or gls, but without success.  Also no luck with
RSiteSearch.

Any suggestions?

I have included below a sample problem for which I think I have the correct
maximum likelihood
estimates for beta and sigma^2 and random effects u., (in res$beta,
res$sig2e and res$u below)
 which I got with a quick hack, but which I would like to reproduce with
some standard R method.

Thanks

Rob

xy <-
structure(c(0.02810951054828, -0.0059751447363, -0.00259486765798,
-0.02408532375275, 0.01233681858059, 0.00335238957434, -0.00437531552244,
0.0085573052387, -0.00295186914992, -5.848326315e-05, 0.05808473441535,
-0.00624215347984, -0.02844949344284, 0.01463634952832, -0.01537450106168,
0.00789132330785, -0.01758284279198, 0.02185240478034, 0.02861913422167,
0.00269505879316, 0.00645232139019, -0.01078894167733, -0.0549443656218,
0.05394083115548, 0.00336404780951, 0.01039785126456, 0.00651310752339,
0.00545016685343, 0.01795693134289, 0.00254833115767, 0.00629453997195,
0.0057251039691, 0.00669739748211, 0.01247983496454, 0.0004798937386,
0.00358346930317, 0.01047475275143, 0.00226774833578, 0.04163824745781,
0.0300726114684, 0.00537941752895, 0.0037613012053, 0.00372696087132,
0.01272218132145, 0.00444575648527, 0.01355607900327, 0.00495450061539,
0.00665815715587, 0.70018484288354, 0.77119868616205, 0.69679575751332,
0.83832303711917, 0.71465870943677, 0.81584990584669, 0.70391345385899,
0.83198374691626, 0.64693819827464, 0.77172704235109, 0.55166943102329,
0.74741880136058, 0.6975613116561, 0.67302720669703, 0.67255775259468,
0.59988447097333, 0.78113776756213, 0.68642481825379, 0.59974397269042,
0.68629725672171, 0.67463529808584, 0.7008583976577, 0.25786272376619,
0.59450936196294, 0.00490892500726, 0.01068001582643, 0.02329564566335,
0.02006917927665, 0.0006186210563, 0.00012864166829, 0.01773028685484,
0.00278397080297, 0.00204658742582, 0.01360133771674, 0.01381488736597,
0.01851231408294, 0.00535066427928, 0.02991886041973, 0.01497999546635,
-0.00826793443986, 0.00515847523844, 0.01793658872856, 0.00964263205073,
0.00625921353626, 0.00390790240252, 0.00237179752867, 0.25453351836804,
0.00663340112391, 0.00010209967539, 8.4110404969007e-06, 4.8349423626399
e-06,
0.00010760136729, 0.00021301765714, 2.203499903e-05, 9.413510752e-05,
4.250192525e-05, 0.00027227682566, 2.106425928e-05, 0.00047273827565,
1.843809247e-05, 1.210348838e-05, 0.00058437528606, 3.159376081803e-06,
5.8875782795912e-06, 0.00016728261028, 0.00010030468982, 0.00049383802885,
0.00025073676102, 3.900598614e-05, 6.284136268e-05, 0.00018436588238,
1.36434877e-05), .Dim = c(24L, 5L), .Dimnames = list(c("1", "2",
"3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14",
"15", "16", "17", "18", "19", "20", "21", "22", "23", "24"),
    c("y", "x1", "x2", "x3", "psi2")))

 xy
               y                         x1                  x2
x3                     psi2
1   2.810951e-02 0.0033640478 0.7001848  0.0049089250 1.020997e-04
2  -5.975145e-03 0.0103978513 0.7711987  0.0106800158 8.411040e-06
3  -2.594868e-03 0.0065131075 0.6967958  0.0232956457 4.834942e-06
4  -2.408532e-02 0.0054501669 0.8383230  0.0200691793 1.076014e-04
5   1.233682e-02 0.0179569313 0.7146587  0.0006186211 2.130177e-04
6   3.352390e-03 0.0025483312 0.8158499  0.0001286417 2.203500e-05
7  -4.375316e-03 0.0062945400 0.7039135  0.0177302869 9.413511e-05
8   8.557305e-03 0.0057251040 0.8319837  0.0027839708 4.250193e-05
9  -2.951869e-03 0.0066973975 0.6469382  0.0020465874 2.722768e-04
10 -5.848326e-05 0.0124798350 0.7717270  0.0136013377 2.106426e-05
11  5.808473e-02 0.0004798937 0.5516694  0.0138148874 4.727383e-04
12 -6.242153e-03 0.0035834693 0.7474188  0.0185123141 1.843809e-05
13 -2.844949e-02 0.0104747528 0.6975613  0.0053506643 1.210349e-05
14  1.463635e-02 0.0022677483 0.6730272  0.0299188604 5.843753e-04
15 -1.537450e-02 0.0416382475 0.6725578  0.0149799955 3.159376e-06
16  7.891323e-03 0.0300726115 0.5998845 -0.0082679344 5.887578e-06
17 -1.758284e-02 0.0053794175 0.7811378  0.0051584752 1.672826e-04
18  2.185240e-02 0.0037613012 0.6864248  0.0179365887 1.003047e-04
19  2.861913e-02 0.0037269609 0.5997440  0.0096426321 4.938380e-04
20  2.695059e-03 0.0127221813 0.6862973  0.0062592135 2.507368e-04
21  6.452321e-03 0.0044457565 0.6746353  0.0039079024 3.900599e-05
22 -1.078894e-02 0.0135560790 0.7008584  0.0023717975 6.284136e-05
23 -5.494437e-02 0.0049545006 0.2578627  0.2545335184 1.843659e-04
24  5.394083e-02 0.0066581572 0.5945094  0.0066334011 1.364349e-05

> res$beta
(Intercept)          x1          x2          x3
  0.1369432  -0.9514077  -0.1672438  -0.5581261
> res$u
            1             2             3             4             5
 5.290481e-03  8.914887e-05 -1.039673e-04 -1.707174e-03  6.828549e-03
            6             7             8             9            10
 6.073888e-04 -2.725684e-03  3.517122e-03 -1.481810e-02  1.257265e-03
           11            12            13            14            15
 1.581453e-02 -4.299002e-04 -2.350822e-03  7.036544e-03  1.467205e-04
           16            17            18            19            20
-1.563975e-04 -7.831411e-03  4.896880e-03  6.726648e-04 -2.295879e-03
           21            22            23            24
-2.078902e-03 -4.358961e-03 -1.027172e-03  1.943738e-03
attr(,"format")
[1] "best"
> res$sig2e
[1] 0.0001721026