Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: rcal, measurement error model question


From   "Austin Nichols" <[email protected]>
To   [email protected]
Subject   Re: st: rcal, measurement error model question
Date   Fri, 12 Dec 2008 11:41:03 -0500

Richard Valliant <[email protected]>:

Below code corrects the aweight typo, and makes better fake variances,
and when I wrote "I am making no use of the
individual-observation-specific error variances you have" I meant
"making no use of the individual-observation-specific error variances
for e2" not e1, since v1 is incorporated in the aweight.

set seed 1
drawnorm e1 v1 e2 v2, n(30) clear
replace v1=exp(v1-2)
replace v2=exp(v2-2)
su e2
loc e2var=r(Var)
su v2
loc v2var=r(mean)
loc r=(`e2var'-`v2var')/`e2var'
eivreg e1 e2 [aw=1/v1], r(e2 `r')
mat r=(`v2var')
rcal (e1=) (w: e2), suuinit(r)

What's the ultimate goal of this exercise?  Is e2 actually supposed to
have a linear effect on e1?

Maybe a simulation would help...


cap prog drop evrcal
prog evrcal, rclass
drawnorm z v1 e2 v2, n(30) clear
replace v1=exp(v1-3)
replace v2=exp(v2-3)
g err1=invnormal(uniform())*v1
g err2=invnormal(uniform())*v2
g o2=e2+err2
g e1=1+(o2)/2+z
g o1=e1+err1
su o2
loc e2var=r(Var)
su v2
loc v2var=r(mean)
loc r=(`e2var'-`v2var')/`e2var'
eivreg o1 o2 [aw=1/v1], r(o2 `r')
return scalar be=_b[o2]
return scalar se=_se[o2]
test o2=.5
return scalar pe=r(p)
mat r=(`v2var')
rcal (o1=) (w: o2), suuinit(r)
return scalar br=_b[w]
return scalar sr=_se[w]
test _b[w]=.5
return scalar pr=r(p)
eret clear
end
simul, reps(10000): evrcal
g rejr=pr<.05
g reje=pe<.05
su, sep(4)



On Fri, Dec 12, 2008 at 10:18 AM, Austin Nichols
<[email protected]> wrote:
> Richard Valliant <[email protected]>:
> I'm pretty sure -rcal- does not do what you want, which I think is
> regressing e1 on e2 while specifying the measurement error variance
> for each observation on both e2 and e1.  It seems from the paper I
> reference that -rcal- will only use the first element of your diagonal
> matrix as the assumed error variance for e2, which you can see
> directly by doing something like:
>
> clear
> drawnorm e1 v1 e2 v2, n(30)
> replace v1=exp(v1/10)
> replace v2=exp(v2/10)
> mkmat v2
> mat D = diag(v2)
> mat d=D[1,1]
> rcal (e1=) (w: e2), suuinit(d)
> rcal (e1=) (w: e2), suuinit(D)
>
> Maybe something like this works for your example (but I hope others on
> Statalist will weigh in):
>
> su e2
> loc e2var=r(Var)
> su v2
> loc v2var=r(mean)
> loc r=(`e2var'-`v2var')/`e2var'
> eivreg e1 e2 [aw=v1], r(e2 `r')
>
> but note that above I am making no use of the
> individual-observation-specific error variances you have.
>
> On Thu, Dec 11, 2008 at 3:03 PM, Austin Nichols <[email protected]> wrote:
>> Richard Valliant <[email protected]>:
>> Maybe you want:
>>  rcal (e1=) (w: e2), suuinit(D)
>> (i.e. leave out only the x not the =x)?
>> See also http://www.stata-journal.com/sjpdf.html?articlenum=st0050
>>
>> On Thu, Dec 11, 2008 at 11:07 AM, Richard Valliant
>> <[email protected]> wrote:
>>> I'm a new user who is trying to fit a simple measurement error model to
>>> a set of estimates from two independent surveys.  The surveys are
>>> measuring the same things.  My data look like
>>>
>>> (e1, v1) = set of 30 estimates and their variances from survey 1
>>> (e2, v2) = set of 30 estimates and their variances from survey 2
>>>
>>> The model I want to fit is basic:
>>> e1 = a + b*e2 + u1
>>> e2 = E2 + u2     (u1 and u2 are the model errors, E2 is E(e2)  )
>>>
>>> I've tried:
>>>  mkmat v2
>>>  mat D = diag(v2)
>>>  rcal (e1) (w: e2), suuinit(D)
>>>
>>> This gives "invalid syntax". If I put some arbitrary variable x in the
>>> model (which I don't want), this works:
>>>    rcal (e1=x) (w: e2), suuinit(D)
>>>
>>> But rcal apparently does not allow aweights to account for v1 =
>>> var(e1).
>>> Is there a way to use rcal or some other procedure to fit the model
>>> above, accounting for the fact that I have (1) estimates of variance for
>>> both e1 and e2 and (2) no covariates measured without error to put in
>>> the model?
>>
>
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2024 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index