Statalist


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

Re: st: Variogram for longitudinal data


From   Austin Nichols <austinnichols@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Variogram for longitudinal data
Date   Tue, 27 Oct 2009 12:29:22 -0400

************************************************
*! xtvariog2 1.1 AN 27 Oct 2009 added panel support
*! based on variog2 1.0.0 NJC 4 March 2005
program xtvariog2, sort
version 8
syntax varlist(numeric min=3 max=3) [if] [in], Width(numlist max=1 >0) ///
 [ Lags(real 0) LIST Generate(string) Id(varlist) noGRaph * ]
if "`id'"=="" {
 xtset
 loc id=r(panelvar)
 }
tempvar ivar t
g long `t'=_n
egen long `ivar'=group(`id')

* check data and options
marksample touse
markout `touse' `ivar'
qui count if `touse'
if r(N) == 0 error 2000
local N = r(N)
if "`generate'" != "" {
 confirm new variable `generate'
 conf new variable `generate'_lag
 }
if `lags' < 0 {
 di as err "lags must be positive"
 exit 198
 }
tempvar ict
egen long `ict'=sum(`touse'), by(`ivar')
su `ict' if `ivar'<., mean
loc maxn=r(max)
local lags = cond(`lags',`lags',int(`maxn'/2))
local lags = min(`lags',`maxn'-5)

* calculation
tempvar lag gamma npairs LAG
tokenize "`varlist'"
args z x y
quietly {
 replace `touse' = -`touse'
 sort `touse' `ivar' `t'
 by `touse' `ivar': gen `lag' = _n if _n<`lags' & `touse'==-1
 label var `lag' "Lag (width `width')"
 gen byte `LAG' = .
 gen double `gamma' = 0
 gen long `npairs' = 0
 label var `gamma' "Semi-variance"
 forval i = 1/`maxn' {
  local I = `i' + 1
  forval j = `I'/`N' {
   by `touse' `ivar':replace ///
    `LAG'=ceil(sqrt((`y'[`i']-`y'[`j'])^2+((`x'[`i']-`x'[`j'])^2))/`width')
   by `touse' `ivar':replace `npairs'=`npairs'+1 if _n==`LAG'
   by `touse' `ivar':replace `gamma'=`gamma'+(`z'[`i']-`z'[`j'])^2 if _n==`LAG'
   }
  }
 replace `gamma' = `gamma' / (2 * `npairs')
 }

* list if desired
if "`list'" == "list" {
 char `lag'[varname] "Lag"
 char `gamma'[varname] "Semi-variance"
 char `npairs'[varname] "# of pairs"
 l `lag' `gamma' `npairs' if !mi(`lag',`gamma'), ///
  subvarname noobs abb(13) sepby(`ivar')
 }

* graph
if "`graph'"=="" {
 local zlab : variable label `z'
 if `"`zlab'"' == "" local zlab "`z'"
 loc o ti(`"Semi-variogram of `zlab'"') i(`ivar') t(`lag') leg(off) `options'
 xtline `gamma' if `gamma'<., ov ysc(r(0 .)) yla(, ang(h)) `o'
 }

* generate if desired
qui if "`generate'" != "" {
 ren `gamma' `generate'
 label var `generate' "Semi-variance of `zlab'"
 ren `lag' `generate'_lag
 }
end
************************************************
use http://www.ats.ucla.edu/stat/stata/faq/ozone
g obs=_n
expand 2
bys station: g year=_n
drop if year==2 & obs>14
set seed 1
replace av8top=av8top+rnormal() if year==2
xtvariog2 av8top lat lon, w(.1) i(year) list g(x)


On Tue, Oct 27, 2009 at 10:00 AM, Stas Kolenikov <skolenik@gmail.com> wrote:
> The standard recommendation applies: take the code of -variog*- and
> modify it to suit your needs (saving a copy under a different name).
> You would have to change all the single observation concepts in
> -variog*- with the panel-level concepts (which is moderately
> difficult), and add another cycle in the computation of `zdsq' scalar
> (which is quite easy).
>
> On 10/26/09, Raphael Fraser <raphael.fraser@gmail.com> wrote:
>> Are there any programs available to estimate variogram for
>>  longitudinal data? Note that -variog- or -variog2- does not apply to
>>  longitudinal data.

*
*   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–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index