I would like Stata to develop some options for programming options for Mata to "easily" develop programs that can use multiple processors at the same time.
-
Login or Register
- Log in with
// Date: 15-Jun-2021
// Use equations shown in Andy Field's document to compute
// Welch's F-test.
// https://www.discoveringstatistics.com/docs/welchf.pdf
// The data AF uses come from this document:
// https://www.discoveringstatistics.com/repository/contrasts.pdf
clear
input byte group y
1 3
1 2
1 1
1 1
1 4
2 5
2 2
2 4
2 2
2 3
3 7
3 4
3 5
3 3
3 6
end
preserve
drop if missing(group)
statsby Ybar=r(mean) Var=r(Var) n=r(N), by(group) clear: summarize y
list
quietly {
generate double Wt = n/Var
generate double WtMean = Wt*Ybar
egen double SumWt = total(Wt)
egen double SumWtMean = total(WtMean)
generate double WelchGM = SumWtMean/SumWt
generate double WtSqDev = Wt*(Ybar-WelchGM)^2
egen double SSmodel = total(WtSqDev)
generate byte k = _N // # of groups
generate double MSmodel = SSmodel/(k-1)
generate double Lterm = (1-Wt/SumWt)^2 / (n-1)
egen double lambda = total(Lterm)
replace lambda = lambda*3/(k^2-1)
generate double MSerror = 1+2*lambda*(k-2)/3
generate double WelchF = MSmodel/MSerror
generate byte df1 = k-1
generate double df2 = 1/lambda
generate double p = Ftail(df1,df2,WelchF)
}
list WelchF-p in 1
display _newline ///
"From SPSS ONEWAY:" _newline ///
"Welch F = "4.320451 _newline ///
" df1 = " 2 _newline ///
" df2 = " 7.943375 _newline ///
" p = " .053738
restore
// Date: 17-Aug-2021
// Compare Welch's F-test to multilevel modeling approach
// in post #3 here:
// https://www.statalist.org/forums/forum/general-stata-discussion/general/1427959-anova-rejection-of-bartlett-test
// Note that I have added the dfmethod(anova) option.
mixed y i.group, ml residuals(independent, by(group)) ///
nolrtest nolog dfmethod(anova)
margins group
margins group, pwcompare(effects) mcompare(noadjust)
// Compare to -regress- with vce(robust)
regress y i.group, vce(robust)
reg log(y) x1 x2 poisson y x x^2 x^3 etc.
reg log(y) x1 x2 poisson y x x^2 x^3 etc.
. sysuse auto
(1978 automobile data)
. expr: reg log(price) weight rep78
Expression _Expr0 := log(price)
-> reg _Expr0 weight rep78
Source | SS df MS Number of obs = 69
-------------+---------------------------------- F(2, 66) = 21.24
Model | 4.00927365 2 2.00463683 Prob > F = 0.0000
Residual | 6.22817243 66 .094366249 R-squared = 0.3916
-------------+---------------------------------- Adj R-squared = 0.3732
Total | 10.2374461 68 .150550678 Root MSE = .30719
------------------------------------------------------------------------------
_Expr0 | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
weight | .000333 .0000513 6.50 0.000 .0002307 .0004354
rep78 | .1272241 .0410658 3.10 0.003 .0452337 .2092146
_cons | 7.196296 .2500145 28.78 0.000 6.697125 7.695466
------------------------------------------------------------------------------
. expr: poisson price weight weight^2 rep78
Expression _Expr0 := weight^2
-> poisson price weight _Expr0 rep78
Iteration 0: log likelihood = -21667.605
Iteration 1: log likelihood = -21655.464
Iteration 2: log likelihood = -21655.463
Poisson regression Number of obs = 69
LR chi2(2) = 36753.75
Prob > chi2 = 0.0000
Log likelihood = -21655.463 Pseudo R2 = 0.4591
------------------------------------------------------------------------------
price | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
weight | -.0005043 .0000145 -34.90 0.000 -.0005326 -.000476
_Expr0 | 1.36e-07 2.20e-09 61.86 0.000 1.32e-07 1.40e-07
rep78 | .0955141 .0018534 51.53 0.000 .0918815 .0991466
_cons | 8.555204 .025696 332.94 0.000 8.50484 8.605567
-
. sysuse auto
(1978 automobile data)
. exprcmd
(expr) . reg log(price) weight rep78
Expression _Expr0 := log(price)
-> reg _Expr0 weight rep78
Source | SS df MS Number of obs = 69
-------------+---------------------------------- F(2, 66) = 21.24
Model | 4.00927365 2 2.00463683 Prob > F = 0.0000
Residual | 6.22817243 66 .094366249 R-squared = 0.3916
-------------+---------------------------------- Adj R-squared = 0.3732
Total | 10.2374461 68 .150550678 Root MSE = .30719
------------------------------------------------------------------------------
_Expr0 | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
weight | .000333 .0000513 6.50 0.000 .0002307 .0004354
rep78 | .1272241 .0410658 3.10 0.003 .0452337 .2092146
_cons | 7.196296 .2500145 28.78 0.000 6.697125 7.695466
------------------------------------------------------------------------------
(expr) . poisson price weight weight^2 rep78
Expression _Expr0 := weight^2
-> poisson price weight _Expr0 rep78
Iteration 0: log likelihood = -21667.605
Iteration 1: log likelihood = -21655.464
Iteration 2: log likelihood = -21655.463
Poisson regression Number of obs = 69
LR chi2(2) = 36753.75
Prob > chi2 = 0.0000
Log likelihood = -21655.463 Pseudo R2 = 0.4591
------------------------------------------------------------------------------
price | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
weight | -.0005043 .0000145 -34.90 0.000 -.0005326 -.000476
_Expr0 | 1.36e-07 2.20e-09 61.86 0.000 1.32e-07 1.40e-07
rep78 | .0955141 .0018534 51.53 0.000 .0918815 .0991466
_cons | 8.555204 .025696 332.94 0.000 8.50484 8.605567
------------------------------------------------------------------------------
(expr) . .
. sysuse auto
(1978 automobile data)
. expr: reg log(price) weight weight^2 rep78
Expression _Expr0 := log(price)
Expression _Expr1 := weight^2
-> reg _Expr0 weight _Expr1 rep78
Source | SS df MS Number of obs = 69
-------------+---------------------------------- F(3, 65) = 18.23
Model | 4.67772488 3 1.55924163 Prob > F = 0.0000
Residual | 5.5597212 65 .085534172 R-squared = 0.4569
-------------+---------------------------------- Adj R-squared = 0.4319
Total | 10.2374461 68 .150550678 Root MSE = .29246
------------------------------------------------------------------------------
_Expr0 | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
weight | -.00067 .0003621 -1.85 0.069 -.0013932 .0000532
_Expr1 | 1.60e-07 5.74e-08 2.80 0.007 4.58e-08 2.75e-07
rep78 | .0963816 .0406237 2.37 0.021 .0152506 .1775127
_cons | 8.768601 .6107286 14.36 0.000 7.548892 9.98831
------------------------------------------------------------------------------
. diest
------------------------------------------------------------------------------
_Expr0 log(price) | Coef. Std. Err. t P>|t|
-----------------------------------------+------------------------------------
weight Weight (lbs.) | -.00067 .0003621 -1.850 0.069
_Expr1 weight^2 | 1.60e-07 5.74e-08 2.796 0.007
rep78 Repair record 1978 | .0963816 .0406237 2.373 0.021
_cons | 8.768601 .6107286 14.358 0.000
------------------------------------------------------------------------------
cap preserve cap drop _all sysuse auto qui sum mpg loc m=r(mean) tw hist mpg, plotr(marg(zero)) xli(`m', lco(red)) saving(gsl1, replace) tw (hist mpg, plotr(marg(zero)) leg(off)) /// (function y=`m', horiz ra(0 .1) lco(red) lpa(solid) saving(gsl2, replace)) drop _all set seed 2345 set obs 2000 gen x=rnormal(0 1) gen y=rnormal(0 1) qui sum x loc mx=r(mean) loc mnx=r(min) loc mxx=r(max) qui sum y loc my=r(mean) loc mny=r(min) loc mxy=r(max) tw scatter y x, plotr(marg(zero)) xli(`mx', lco(red)) yli(`my',lco(red)) xla(-2(1)5) saving(gsl3, replace) tw (scatter y x, plotr(marg(zero)) xsc(r(`mnx' `mxx')) xla(-2(1)5) ysc(r(`mny' `mxy')) leg(off)) /// (function y=`my', ra(`mnx' `mxx') lco(red) lpa(solid) ) /// (function y=`mx', horiz ra(`mny' `mxy') lpa(solid) lco(red) saving(gsl4, replace)) cap restore
Comment