我的同事Liwen

Liwen是我们公司marketing部门strategic planning team的人。在我的印象中，她是一个非常关注自己职业发展和周围环境的人。记得她上班的第一天就逐个到每个人的cub里面进行自我介绍，还问我们是不是Excel的大牛，以后要向我们请教。不仅如此，她在networking上面非常积极，刚上班的第三个星期就给我们发信要组织marketing跟distribution的中国人lunch and learn，而且每次lunch and learn之前她都会精心准备一番，画小map给第一次来我们会议室的人，Lunch and learn之后的会议总结等等。在我印象中，她每个星期都要跟几个不同的人吃饭，在很短的时间内，她就认识了不少人。 她特别关心我们周围的人都在做什么，经常问的一句话就是，“你在做什么project?”或者“你跟你老板处的怎么样啊？”。正因为如此，大家都比较喜欢跟她交流自己的工作，发发牢骚，聊聊天什么的。她会很仔细观察周围的人，然后跟我说，谁谁谁比较会做事情啊，感觉谁谁谁的脾气很好啊，技术很好啊之类的。

店靠衣装

今天跟好朋友QQ去H-mark，是一个韩国超市。这个超市给我的明显的一个感觉是干净，整洁。整个超市的面积跟Super 88差不多，但是明显能感觉到各个section安排的极有条理，每个section还有一些店员做展示，我们就是因为看了店员展示的食物买了两大盒牛肉。韩国不是有非常多的小菜吗？超市非常精心的把每个小菜打包，封装好，甚至连肉类也包装得给你感觉是打开来就能吃，而且美美的透明的包装（有的包装袋上面还点缀点小绿叶）带带给人一种感觉：我们的菜非常新鲜哦～。相比之下，super 88，东西摆放不仅凌乱，食物表面也经常油腻腻，让人想拿起来仔细看看的欲望都没有。甚至还经常会买到过期食物。

Super 88地点在Boston的黄金地带，客流量非常大，开车只要15分钟。而H-mart在非常偏远的郊区，开车过去都要40分钟之久。但是H-mart却是我经常去逛并且非常喜欢逛的一家店。其实super 88本有更好的条件能有更多的客户，却因为经营不用心，很多客户都宁愿到别的小店购买更新鲜的蔬菜。

做为一家每天都为消费者提供日常用品，消费食物的店，客户一般比较在意的是：(1) 蔬菜水果生鲜肉类的新鲜程度 (2) 价格 (3) 购物环境 (4) 产品外包装 (5) check out款台的服务。

今天跟表弟打电话，他还在找工作阶段，说现在面试的机会比较少了。找工作也是分周期的，一般到年底机会就是比平时要少一点。美国也是这样的。为了能够控制他打游戏的时间，我跟他约定下周四或者周五提交一份他的作品（他是做动漫设计的哦，很酷的职业吧），而且要比上次给我的那份还要好，期待ing!

最近在读：《自控力》作者Kelly McGonigal, Ph.D

提高自控能力的第一步：学会每天做5-15分钟冥想（刚开始可以是5分钟，后来尽快拓展到10-15分钟）。

1.  原地不动，安静坐好。坐在椅子上，双腿平放在地上或者盘坐在垫子上。背挺直，双腿放在膝盖上，冥想时一定不能烦躁，这是自控力的基本保证。
2. 注意你的呼吸。吸气时在脑海中默念“吸”，呼气时在脑海中默念“呼”。当你发现自己走神的时候，将注意力重新集中在呼吸上。这种反复的训练，能够让大脑前皮质处在高速开启的模式，让大脑处理压力和冲动的区域更加稳定。
3.   感受呼吸，弄清楚自己是怎么走神的。几分钟后你可以不再专注于默念“呼”或者“吸”而是专注于呼吸本身，你会注意到空气从鼻子和嘴巴进入和呼出的感觉，感觉到吸气时胸腹部的扩张跟呼气时胸腹部的收缩。

   冥想不是让你什么都不想，而是让你不要太分心，不要忘记了最初的目标，如果你在冥想的时候无法集中注意力，别担心，你只需多做练习，将注意力重新集中在呼吸上。  
  核心思想：意志力实际上是“我要做”，“我不要”，“我想做”三种力量，他们协同努力，让我们变成更好的自己。

I'm Back!

居然已经有2年没有写过博客了，用“居然”这个词也不准确，其实自己也知道。原因嘛，各种：主要是懒。懒得去想太多，更别提写了，想总是比写快吧。有时候事情在心里一过就过去了，回顾，大抵也是件费尽的事情。写，要落笔，有时候竟然不知用什么词比较好(很久没写中文了)，还好打字方式有拼音联想，要不真是错别字一大堆(估计这就是为啥现在那个“汉字英雄”，“汉字拼写大会”这类节目那么火了)。

好了，不解释鸟，装得像大明星有一帮粉丝望眼欲穿的等着看你博客更新似的。老娘想写了，不行么？

Parameter Expansion in Bayesian Hierarchical Modeling

This is the project that I did in my Bayeisan Statistics course, below are the snapshots of my part of my paper.

Abstract: Hierarchical model is devoted to facilitate the simultaneous estimation of severalparameters over similar units. However, some problems pertinent to Bayesian hierarchicalmodeling remain unsolved, that is: if the standard deviation for the second layer of hierarchicalmodel (also called between-study standard deviation) has broad peak at zero, somenoninformative prior such as flat uniform prior, IG(0.01,0.01) which are normally adopted inresearch, may lead to insensitivity in the estimation of such smoothing variance. Apart fromthis, in this prior setting, convergence based on EM and Gibbs sampling may become lessconvergence. In this paper, we bring forward a multiplicative parameter-expansion methodto reparameterize hierarchical model in the context of Bayesian inference which facilitatesconvergence and possesses decent properties. Illustrations in terms of simulation will bedelivered to reveal the two-fold essences of this method.

Require for paper"Parameter Expansion in Bayeisan Hierarchical Modeling" ? Simply send an email to stefanie.cao@gmail.com with title "require paper hierarchical modeling".

Bayesian Changing-Point Analysis on US Stock Price During Financial Crisis

This is the project that I have done in my Bayesian Statistics Course, below are snapshots of my paper.

Abstract: Breaks in stock market are usually motivated by an exogenous changein surrounding economy uctuation that precipitates a change in regres-sion regimes. However, di erent industries react di erently to economyfactors (such as unemployment,ination,prime interest and oil price)bothin response time and degree of severity. In this article, Bayesian Changing-Point analysis has been used to detect change point in stock price in fourindustries ranging from automobile, nance, hi-tech and fast moving con-sume goods. The time span is set to be Jan.1,2007 to Dec.30,2008 whichis commonly believed cover the whole process of nancial crisis. With thismodel, we can check the change point location and corresponding poste-rior probability. Inferences for the regression coeffcients before and afterchanging points indicate alternation in sensitivity to economy factors inthose industries.

Require report: send me email via stefanie.cao@gmail.com with title "Bayesian Change Point"

Is p value telling the truth?

Consider a series of experimental testing in drug efficiency, denoted by D1,D2, ...,D20.Suppose now we have two hypotheses:H0 : Drug is not efficient. H1 : Drug is efficient.if cut-off point is preset as 0.05, suppose one of the test results is a p.value=0.045, andanother is p.value=0.016, we statistically reject H0, drawing conclusion that both drugsare efficient. As mentioned above, we need to do 20 tests under the same hypothesisstructure. Then, basically, we can obtain 20 p values, suppose all the p values can belisted below:
Table 1 Drug TestingDrug
D1 D2 D3 D4 D5 D6 D7 D8 D9 D10
p value 0.41 0.31 0.049 0.045 0.016 0.21 0.30 0.209 0.102 0.122
Drug D11 D12 D13 D14 D15 D16 D17 D18 D19 D20
p value 0.121 0.003 0.091 0.40 0.273 0.192 0.167 0.311 0.28 0.22
But the problem is how strong is the evidence that the non-efficient drug istruly coming from non-efficient group?
For more details, see my report, require it ? just send email to stefanie.cao@gmail.com with title "require for p value report".

Is the dataset coming from binomial trials? (let bayesian check)

Source: Bayesian Data Analysis, Andrew Gelman, John B.Carlin etc. (2004, second edition), page 163

Alright, let's start our business.
Consider two dataset: dataset A =(1,1,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0), and dataset B=(1,0,1,0,0,0,1,1,0,0,1,1,1,1,0,0,1,0,1,1).
Question: given dataset A, can we make inference on "theta", the probability that a flipping coin ends up with face.
How does bayesian execute inference?
step 1: assume dataset A comes from binomial distribution (look, this is usually given in textbook problem, in reality, it is not necessarily true, however, people tend to take for granded that dataset A follows binomial distribution)
step 2: figure a prior for "theta". maybe uniform in most cases. make sense, right?
step 3: with prior and likelihood, bayesian can deduce the posterior distribution of "theta".
step 4: done. Give me five....!

However, the dataset A has been modeled as a specified number of iid Bernoulli trials with a uniform prior distribution on the probability of success, say, theta may not actually follow preassumed distribution. I did it a lot without giving even one second of thinking potential assumption really holds. The observed autocorrelation on dataset A is evidence that the model is flawed. To quantify the evidence, we can perform a posterior predictive distributino of T(y^{rep}) by simulation, that is, we assume that dataset A follows iid Bernoulli trials, then we calculate the posterior distribution for the switch.number (the switch.number is the number of times that data change from 0 to 1, or 1 to 0. either way.) and then caculate the probability that posterior switch.number greater than the actually switch.number (in the dataset A , the switch.number is 3). Let's do the experiment! (just for dataset A, as for dataset B, readers, you are smart enough to DIY!)

Following is histogram of simulation of posterior distribution of switch.number.

(Graphic is too small to see ?? just click it, the enlarged one will show in another window)
R program:
y <- c(1,1,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0) # this is the orginal data which is assumed to be binomial distributed, however, whether assumption is valid still pending, leading us to investiage via following program:

set.seed(40)
theta.post <- rbeta(1,sum(y)+1,length(y)-sum(y)+1)

set.seed(40)

y.rep <- array(rbinom(2000,1,theta.post),c(20,100))

switch.number <- array(0,100,1)

for (j in 1:100)

{

for (i in 1:19)

{

if (y.rep[i,j]!=y.rep[i+1,j])

switch.number[j] <- switch.number[j]+1}

}

hist(switch.number,probability=TRUE,breaks=c(1.5:13.5),main = "posterior predictive distribution of number of switches",xlab="switch.number")

lines(density(switch.number),col="red",lwd=3)

abline(v=3, lwd=3,col="blue")

p.value <- sum(switch.number<=3)/100

p.value

result: p.value=0.02 indicating that there is no adequate evidence support the null hypothesis that dataset A comes from binomial distribution.

R function to generate random multivariate Gaussians

# a function to generate random multivariate Gaussians
rmultnorm <- function(n, mu, vmat, tol = 1e-07)
{
p <- ncol(vmat)
if(length(mu)!=p)
stop(paste("mu vector is the wrong length:",length(mu)))
if(max(abs(vmat - t(vmat))) > tol)
stop("vmat not symmetric")
vs <- svd(vmat)
vsqrt <- t(vs$v %*% (t(vs$u) * sqrt(vs$d)))
ans <- matrix(rnorm(n * p), nrow = n) %*% vsqrt
ans <- sweep(ans, 2, mu, "+")
dimnames(ans) <- list(NULL, dimnames(vmat)[[2]])
ans
}

Variance Reduction Using the Antithetic Variables Scheme

Most monte carlo estimates are based on iid samples, but algorithms based on negatively correlated samples maybe superior in terms of producing estimates with smaller variance.Suppose, we use monte carlo method to approximate mean, second moment, median and 75% quantile of a distribution with corresponding denstiy function: 2*X^3 I(x>1). The inverse method is usually used to generate samples follows such distribution. However, if we use both 1-U1 and U1 instead of U1. Here is R program to present variance reduction using both 1-U1 and U1.f(x)=2*X^3I(x>1)

(unfinish...)

R code to solve Cubic Equation

This is create function to solve cubic equation.

note: default form of a cubic equation in this "cubicEq" function is: a*X^3+b*X^2+c*X+d=0. 

cubicEq<-function(a,b,c,d)

{  

   x1<-c()

   x2<-c()

   x3<-c()

   for ( i in 1:2)

   {

      f<-(3*c[i]/a[i]-(b[i]^2/(a[i]^2)))/3

      g<-(2*b[i]^3/a[i]^3-(9*b[i]*c[i]/a[i]^2)+27*d[i]/a[i])/27

      h<-(g^2/4)+(f^3/27)

      if (h<=0)

         {

             count<-(g*g/4-h)^(1/2)

             j<-count^(1/3)

             k<-acos(-g/(2*count))

             L<-j*(-1)

             M<-cos(k/3)

             N<-3^(1/2)*sin(k/3)

             p<-(-1)*(b[i]/(3*a[i]))

             x1[i]<-2*j*cos(k/3)-b[i]/(3*a[i])

             x2[i]<-L*(M+N)+p

             x3[i]<-L*(M-N)+p

          }

       if ( h>0 )

          {

             R<- -g/2+h^(1/2)

             S<-R^(1/3)

             T<- -g/2-h^(1/2)

             U<--(-T)^(1/3)

             x1[i]<-(S+U)-(b[i]/(3*a[i]))

             x2[i]<- complex(real=-(S+U)/2-b/(3*a[i]),imaginary=((S-U)*sqrt(3)/2))

             x3[i]<- complex(real=-(S+U)/2-b/(3*a[i]),imaginary=-((S-U)*sqrt(3)/2))

          }

      }

return(list(x1=x1,x2=x2,x3=x3))

}

R code to generate symmetric matrix M

Note: this is a minic program to generate a symmetric matrix "M".
nRow-----the number of rows in M
nCol------the number of columns in M
diag-------diagonal elements of M
nondiag----non-diagonal elements of M

matr<-function(nRow,nCol,diag,nondiag)
{M<-array(,dim=c(nRow,nCol))
for ( i in 1:nRow)
{ for ( j in 1:nCol)
{
if (i!=j)
M[i,j]<-nondiag else M[i,j]<-diag
j+1
}
i+1
}
return(M)
}

refuse messy

the less stuff on your desk, the higher efficiency you will get.

Experiment Design on Aluminum Alloy Research

This is a mini project from my DOE course.
Problem Description:
A structural engineer is studying the strength of aluminum alloy purchased from three vendors. Each vendor submits the alloy in standard-sized bars of 1.0, 1.5 and 2.0 inches. The processing of different sizes of bar stock from a common ingot involves different forging techniques, and so this factor may be important. Furthermore, the bar stock is forged from ingots made in different heats. Each vendor submits two test specimens of each size bar stock from three heats. The resulting strength data follow. Analyze the data,assuming that vendors and bar size are fixed and heats are random. Use the resticted model and unrestricted model of the mixed model.

Data structure:


report snapshot:
Report Requet: send email to stefanie.cao@gmail.com with title "mixed model report request", you will receive report within three days.

Exploring Female Body Weight using Multiple Linear Regression

This is a mini project from my linear model course.
1. Introduction
Human Body weight varies according to different human group. In this project; we are interested in human body weight in female within the age of 27. A weight equation for the female group in terms of body builds variables was found using multiple regression technique. Also series of T test have been executed to verify the relative significant importance of each variable, determining variable entries. The whole procedure can be described as follow figure1-1.
2. Procedure

3. Data

The initial measurements were taken at San Jose State University and at the U.S. Naval Postgraduate School in Monterey, California. These were primarily female individuals with in 27 years old, all physically active (several hours of exercise a week). Nine skeletal measurements (diameter measurements) were included in the study. Twelve girth (or circumference) measurements are used in the study. They, in contrast to skeletal (diameter) measurements, are not fixed over time except for the three "bony" girths of the wrist, knee, and ankle, which remain relatively constant over the life span. The other nine girths, the changeable ones, were measured at these sites: shoulder, chest, waist, navel, hip, thigh, bicep, forearm, and calf. In addition to the skeletal and girth measurements, recorded in centimeters (cm) is height (cm). Data information stem from such website: http://www.sci.usq.edu.au/staff/dunn/Datasets /applications/biology/body.html

4. Result

...

5. Conclusion

...

6. Appendix

code:

PROC IMPORT OUT= SASUSER.BODY
DATAFILE= "C:\Documents and Settings\Administrator\desktop\Project\body.dat";
DBMS=TAB REPLACE;
GETNAMES=YES;
DATAROW=2;
RUN;

/*----------------------------------------

sifting data: gender=female, age less than 27

-----------------------------------------*/
proc iml;
use sasuser.body ;
read all var {'biacromial' 'pelvic_breadth' 'bitrochanteric' 'chest_depth' 'chest_diam' 'elbow_diam' 'wrist_diam'
'knee_diam' 'ankle_diam' 'shoulder_girth' 'chest_girth' 'waist_girth' 'navel_girth' 'hip_girth' 'thigh_girth' 'bicep_girth' 'forearm_girth' 'knee_girth' 'calf_girth' 'ankle_girth' 'wrist_girth' 'height' 'age' }into x;
read all var {'weight'} into y;
run;

data sasuser.bodysl;
set sasuser.body;
if 0 le age le 27 & gender='1'; run;

/*----------------------------
Normality Test of Weigth
-----------------------------*/

proc univariate data=sasuser.bodysl normal;
var weight;
histogram/normal;
run;

/*---------------------------------
Correlation between variables
----------------------------------*/

proc corr data=sasuser.bodysl;
var biacromial pelvic_breadth bitrochanteric chest_depth chest_diam elbow_diam wrist_diam knee_diam ankle_diam shoulder_girth chest_girth waist_girth navel_girth hip_girth thigh_girth
bicep_girth forearm_girth knee_girth calf_girth ankle_girth wrist_girth height;
run;

data sasuser.body;
set sasuser.bodysl;
keep weight biacromial pelvic_breadth bitrochanteric chest_depth wrist_diam knee_diam ankle_diam chest_girth thigh_girth calf_girth ankle_girth height;
run;

/*----------------------------------
Initial multiple regression Model
-----------------------------------*/
proc iml;
use sasuser.bodysl;
read all var {'biacromial' 'pelvic_breadth' 'bitrochanteric' 'chest_depth' 'wrist_diam' 'knee_diam' 'ankle_diam' 'chest_girth' 'thigh_girth' 'calf_girth' 'ankle_girth' 'height' } into x;
read all var{'weight'} into y;
n=nrow(x);
m=ncol(x);
x=j(n,1,1)x;
beta=inv(x`*x)*x`*y;
print beta;
run;

/*-------------------------------------
First round t-test for coefficient
--------------------------------------*/
proc iml;
use sasuser.bodysl;
read all var {'biacromial' 'pelvic_breadth' 'bitrochanteric' 'chest_depth' 'wrist_diam' 'knee_diam' 'ankle_diam' 'chest_girth' 'thigh_girth' 'calf_girth' 'ankle_girth' 'height' } into x;
read all var{'weight'} into y;
n=nrow(x);
m=ncol(x);
x=j(n,1,1)x;
beta=inv(x`*x)*x`*y;
print beta;
sse=y`*y-beta`*x`*y;
df_sse=n-m-1;
mse=sse/df_sse;
varbeta=inv(x`*x)*mse;
print varbeta;
do i=2 to 13;
g=varbeta[i,i];
sqr=sqrt(g);
ttest1=beta[i,1]/sqr;
if ttest1<0 p="2*probt(ttest1,89);">0 then p=2*(1-probt(ttest1,89));
print p;
end;
run;

/*------------------------------------------------------------------
P value of variable Pelvic_breadth=0.7289, drop pelvic_breadth
------------------------------------------------------------------*/
proc iml;
use sasuser.bodysl;
read all var {'biacromial' 'bitrochanteric' 'chest_depth' 'wrist_diam' 'knee_diam' 'ankle_diam' 'chest_girth' 'thigh_girth' 'calf_girth' 'ankle_girth' 'height' } into x;
read all var{'weight'} into y;
n=nrow(x);
m=ncol(x);
x=j(n,1,1)x;
beta=inv(x`*x)*x`*y;
sse=y`*y-beta`*x`*y;
df_sse=n-m-1;
mse=sse/df_sse;
varbeta=inv(x`*x)*mse;
print varbeta;
do i=2 to 12;
g=varbeta[i,i];
sqr=sqrt(g);
ttest1=beta[i,1]/sqr;
if ttest1<0 p="2*probt(ttest1,90);">0 then p=2*(1-probt(ttest1,90));
print p;
end;
run;

/*-----------------------------------------------------------
P value of variable Wrist_diam=0.49385, drop Wrist_diam
-----------------------------------------------------------*/
proc iml;
use sasuser.bodysl;
read all var {'biacromial' 'bitrochanteric' 'chest_depth' 'knee_diam' 'ankle_diam' 'chest_girth' 'thigh_girth'
'calf_girth' 'ankle_girth' 'height' } into x;
read all var{'weight'} into y;
n=nrow(x);
m=ncol(x);
x=j(n,1,1)x;
beta=inv(x`*x)*x`*y;
sse=y`*y-beta`*x`*y;
df_sse=n-m-1;
mse=sse/df_sse;
varbeta=inv(x`*x)*mse;
print varbeta;
do i=2 to 11;
g=varbeta[i,i];
sqr=sqrt(g);
ttest1=beta[i,1]/sqr;
if ttest1<0 p="2*probt(ttest1,91);">0 then p=2*(1-probt(ttest1,91));
print p;
end;
run;

/*------------------------------------------------
P value of variable Height=0.58, drop height
-------------------------------------------------*/
proc iml;
use sasuser.bodysl;
read all var {'biacromial' 'bitrochanteric' 'chest_depth' 'knee_diam' 'ankle_diam' 'chest_girth' 'thigh_girth' 'calf_girth' 'ankle_girth'} into x;
read all var{'weight'} into y;
n=nrow(x);
m=ncol(x);
x=j(n,1,1)x;
beta=inv(x`*x)*x`*y;
sse=y`*y-beta`*x`*y;
df_sse=n-m-1;
mse=sse/df_sse;
varbeta=inv(x`*x)*mse;
print varbeta;
do i=2 to 10;
g=varbeta[i,i];
sqr=sqrt(g);
ttest1=beta[i,1]/sqr;
if ttest1<0 p="2*probt(ttest1,92);">0 then p=2*(1-probt(ttest1,92));
print p;
end;
run;

/*-------------------------------------------------------------
P value of variable ankle_girth=0.6952, drop ankle_girth
-------------------------------------------------------------*/
proc iml;
use sasuser.bodysl;
read all var {'biacromial' 'bitrochanteric' 'chest_depth' 'knee_diam' 'ankle_diam' 'chest_girth' 'thigh_girth' 'calf_girth' } into x;
read all var{'weight'} into y;
n=nrow(x);
m=ncol(x);
x=j(n,1,1)x;
beta=inv(x`*x)*x`*y;
sse=y`*y-beta`*x`*y;
df_sse=n-m-1;
mse=sse/df_sse;
varbeta=inv(x`*x)*mse;
print varbeta;
do i=2 to 9;
g=varbeta[i,i];
sqr=sqrt(g);
ttest1=beta[i,1]/sqr;
if ttest1<0 p="2*probt(ttest1,93);">0 then p=2*(1-probt(ttest1,93));
print p;
end;
run;
/*----------------------
Confidence Interval
-----------------------*/
proc iml;
use sasuser.bodysl;
read all var{'biacromial' 'bitrochanteric' 'chest_depth' 'knee_diam' 'ankle_diam' 'chest_girth' 'thigh_girth'
'calf_girth' } into x;
read all var{'weight'} into y;
t=tinv(0.975,93);
print t;
n=nrow(x);
m=ncol(x);
x=j(n,1,1)x;
beta=inv(x`*x)*x`*y;
sse=y`*y-beta`*x`*y;
df_sse=n-m-1;
mse=sse/df_sse;
varbeta=inv(x`*x)*mse;
upperb=j(9,1);
lowerb=j(9,1);
do i=1 to 9;
g=varbeta[i,i];
sqr=sqrt(g);
upperb[i,1]=beta[i,1]+t*sqr;
lowerb[i,1]=beta[i,1]-t*sqr;
end;
print upperb;
print lowerb;
run;

/*------------------------------
ANOVA Table Using Proc IML
-------------------------------*/
proc iml;
use sasuser.bodysl;
read all var {'biacromial' 'bitrochanteric' 'chest_depth' 'knee_diam' 'ankle_diam' 'chest_girth' 'thigh_girth' 'calf_girth' } into x;
read all var{'weight'} into y;
n=nrow(x);
m=ncol(x);
vector=j(102,8);
do i=1 to 8;
xbar=mean(x[:,i]);
do count=1 to 102;
vector[count,i]=xbar;
end;
end;
print vector;
cenX=j(102,8);
do k=1 to 102;
do q=1 to 8;
cenX[k,q]=x[k,q]-vector[k,q];
end;
end;
print cenX;
ybar=mean(y[:,1]);
print ybar;
cenbeta=inv(cenX`*cenX)*cenX`*y;
ssr=cenbeta`*cenX`*y;
sse=y`*y-n*(ybar*ybar)-cenbeta`*cenX`*y;
sst=y`*y-n*ybar*ybar;
df_sse=n-m-1;
df_ssr=m;
mse=sse/df_sse;
print ssr;
msr=ssr/df_ssr;
f=msr/mse;
print beta;
print sse ssr df_sse df_ssr;
print mse msr;
print f;
pvalue=1-probf(125.008,8,93);
print pvalue;
run;

/*----------------------------
ANOVA Table Using PROC REG
----------------------------*/
proc reg data=sasuser.body ;
model weight=biacromial bitrochanteric chest_depth knee_diam ankle_diam chest_girth thigh_girth calf_girth/clb;
run;

Request for final report, please send email to stefanie.cao@gmail.com with title: ''Request for Human Weight Research Report''.

PROC ANOVA and PROC GLM Summary

This is meant to be a brief summary of the syntax of the most widely used statements with PROC ANOVA and PROC GLM. There are actually more statements and options that can be used with proc ANOVA and GLM -- you can find out by typing HELP GLM in the command area on the main SAS Display Manager Window.
In the statements below, uppercase is used for keywords, lowercase for things you fill in. Variable names are no more than 8 chars. in length.
PROC ANOVA
PROC ANOVA handles only balanced ANOVA designs. PROC ANOVA DATA=datasetname;
CLASS factorvars;

MODEL responsevar = factorvars; /* See below */
MEANS factorvars / BON /* Bonferroni t-tests,
g=r(r-1)/2 */
T or LSD /* Unprotected t-tests */
TUKEY /* Tukey studentized range */
SCHEFFE /* Scheffe contrasts */
ALPHA=pvalue /* default: 5% */
CLDIFF /* Confidence limits */
LINES /* Non-significant subsets */
; /* use LINES or CLDIFF */

PROC GLM
PROC GLM handles any ANOVA/regression/ANCOVA design. PROC GLM DATA=datasetname;
CLASS factorvars;
e.g. CLASS A B SEX;
MODEL responsevar = factorvars
/ options ; /* Not needed yet */
RANDOM factorvars / TEST; /* If any random factors, list
then here (after MODEL) */
TEST H=effects E=effect; /* To specify an error term other
than the residual MS */
eg, TEST H=A B E=AB; /* 2-way design with A,B random */
MEANS factorvars / options ; /* same as for ANOVA */
LSMEANS factorvars /* Least squares & adjusted
means for ANCOVA */
/ STDERR /* .. and std errors */
PDIFF ; /* ... and p-values for diff */
CONTRAST 'label' factor weights ;
eg, CONTRAST 'Linear' SUGAR -3 -1 1 3 ;
CONTRAST 'Quad ' SUGAR 1 -1 -1 1 ;
ESTIMATE 'name' effect values... / options; /* Only with GLM */
The ESTIMATE statement constructs and tests linear combinations
(predicted values and contrasts) of the parameters.
eg, ESTIMATE 'A1 vs A2' A 1 -1 0 0 / divisor=2;
ESTIMATE 'A2 vs A3,4' A 0 2 -1 -1 / divisor=2;
OUTPUT OUT=datasetname P=fitvar /* Predicted values */
R=residvar ; /* Residuals */
REPEATED factorname levels(levelvalue) contrast;
Models
These illustrate types of MODEL statements that ANOVA and GLM can handle. (Assume A, B, C are class variables; X1, X2, X3 are quantitative, regression variables) MODEL Y = X1; /* Simple linear regression */
MODEL Y = X1 X2 X3; /* Multiple regression */
MODEL Y = X1 X1*X1 X1*X1*X1; /* Polynomial regression */
MODEL Y = A; /* One way anova */
MODEL Y = A B; /* Two-way, main effects only */
MODEL Y = A B A*B; /* Two-way, factorial with
interaction */
MODEL Y = A B; /* Two-way, same as above */
MODEL Y = A B C A*B A*C /* Three-way, complete */
B*C A*B*C; /* factorial */
MODEL Y = A B C; /* The same, using "" notation */
Contrasts
The following table gives coefficients for contrasts to represent linear, quadratic, etc. trend of a quantitative factor on a CONTRAST statement. E.g., for a 3-level factor, use CONTRAST 'linear' DELAY -1 0 1;
CONTRAST 'quad' DELAY 1 -2 1;
Contrast coefficients for Trend Analysis
(Valid when X levels are equally spaced & N's are equal) Coefficients, c(i)
r Trend X=1 2 3 4 5 6 7 sum c(i)**2
---------------------------------------------------------
3 Linear -1 0 1 2
Quad 1 -2 1 6
--------------------------------------------------
4 Linear -3 -1 1 3 20
Quad 1 -1 -1 1 4
Cubic -1 3 -3 1 20
--------------------------------------------------
5 Linear -2 -1 0 1 2 10
Quad 2 -1 -2 -1 2 14
Cubic -1 2 0 -2 1 10
Quartic 1 -4 6 -4 1 70
--------------------------------------------------
6 Linear -5 -3 -1 1 3 5 70
Quad 5 -1 -4 -4 1 5 84
Cubic -5 7 4 -4 -7 5 180
Quartic 1 -3 2 2 -3 1 28
--------------------------------------------------
7 Linear -3 -2 -1 0 1 2 3 28
Quad 5 0 -3 -4 -3 0 5 84
Cubic -1 1 1 0 -1 -1 1 6
Quartic 3 -7 1 6 1 -7 3 154
--------------------------------------------------
8 Linear -7 -5 -3 -1 1 3 5 7 168
Quad 7 1 -3 -5 -5 -3 1 7 168
Cubic -7 5 7 3 -3 -7 -5 7 264
Quartic 7 -13 -3 9 9 -3 -13 7 616
--------------------------------------------------
9 Linear -4 -3 -2 -1 0 1 2 3 4 60
Quad 28 7 -8 -17 -20 -17 -8 7 28 2772
Cubic -14 7 13 9 0 -9 -13 -7 14 990
Quartic 14 -21 -11 9 18 9 -11 -21 14 2002
--------------------------------------------------

Project2: (County Demographic Information)CDI research

This is my test project, though incomplete, still can be your reference for study and research.
Introduction:
The general purpose of Multiple Regressions (the term was first used by Pearson, 1908) is to learn more about the relationship between several independent or predictor variables and a dependent or criterion variable. In practical application, variables can be numerous, which make it inconvenient to build regression model by all variables (may introduce Multi-Colinearity), therefore, alternatives models which use just some variables out of the whole variables are usually brought forward, an optimal model will be chosen from alternative models. In project 1, two alternative models for predicting the number of active physicians in county demographic information (CDI) have been built and several diagnoses been made to compared the model’s fitness. Besides that, by building Multiple Regression Models according to different sample population, we can study the relations between different sample populations and dependent variable, exploring further how sample population influence variance on dependent variable. Project 2 is a such kind of case that concerns relations between the numbers of serious crimes in CDI against three predictors in four regions (four difference sample populations).
data source:
http://apps.csom.umn.edu/nachtsheim/Kutner/Appendix%20C%20Data%20Sets/APPENC02.txt
snapshot of my report:

If you intrested in CDI Research Report, please send a request email to stefanie.cao@gmail.com with email title of "Request for CDI report". I promise to send you report within in three days.

Project: Credit Scoring Via Logistic Model

Introduction:
Logistic regression model is often applied to large finacial databases. For example, credit scoring is a method of modeling the influence of predictors on the probability that a consumer is credit worthy. The data archive found is in: http://www.stat.uni-muenchen.de/service/datenarchiv/kredit/kredit.asc.
Take a snapshot of my reprot:

Structure of my report:
Abstract.............................................................1
section 1: project objective......................................................................4
section 2: Project Procedure.................................................................. 4
section 3: Selection Algorithms and Variable Description..................5
sectoin 4: Reference Model Building and Diagnose.............................7
section 5: Final Model Justification......................................................14
section 6: conclusion...............................................................................15
supplement.........................................15

Want to get final study report?
Send email to stefanie.cao@gmail.com with title of report request on credit scoring. promise you get the report within three days.

What is the difference between categorical, ordinal and interval variables?

In talking about variables, sometimes you hear variables being described as categorical (or sometimes nominal), or ordinal, or interval. Below we will define these terms and explain why they are important.

Categorical

A categorical variable (sometimes called a nominal variable) is one that has two or more categories, but there is no intrinsic ordering to the categories. For example, gender is a categorical variable having two categories (male and female) and there is no intrinsic ordering to the categories. Hair color is also a categorical variable having a number of categories (blonde, brown, brunette, red, etc.) and again, there is no agreed way to order these from highest to lowest. A purely categorical variable is one that simply allows you to assign categories but you cannot clearly order the variables. If the variable has a clear ordering, then that variable would be an ordinal variable, as described below.

Ordinal

An ordinal variable is similar to a categorical variable. The difference between the two is that there is a clear ordering of the variables. For example, suppose you have a variable, economic status, with three categories (low, medium and high). In addition to being able to classify people into these three categories, you can order the categories as low, medium and high. Now consider a variable like educational experience (with values such as elementary school graduate, high school graduate, some college and college graduate). These also can be ordered as elementary school, high school, some college, and college graduate. Even though we can order these from lowest to highest, the spacing between the values may not be the same across the levels of the variables. Say we assign scores 1, 2, 3 and 4 to these four levels of educational experience and we compare the difference in education between categories one and two with the difference in educational experience between categories two and three, or the difference between categories three and four. The difference between categories one and two (elementary and high school) is probably much bigger than the difference between categories two and three (high school and some college). In this example, we can order the people in level of educational experience but the size of the difference between categories is inconsistent (because the spacing between categories one and two is bigger than categories two and three). If these categories were equally spaced, then the variable would be an interval variable.

Interval

An interval variable is similar to an ordinal variable, except that the intervals between the values of the interval variable are equally spaced. For example, suppose you have a variable such as annual income that is measured in dollars, and we have three people who make $10,000, $15,000 and $20,000. The second person makes $5,000 more than the first person and $5,000 less than the third person, and the size of these intervals is the same. If there were two other people who make $90,000 and $95,000, the size of that interval between these two people is also the same ($5,000).

Why does it matter whether a variable is categorical, ordinal or interval?

Statistical computations and analyses assume that the variables have a specific levels of measurement. For example, it would not make sense to compute an average hair color. An average of a categorical variable does not make much sense because there is no intrinsic ordering of the levels of the categories. Moreover, if you tried to compute the average of educational experience as defined in the ordinal section above, you would also obtain a nonsensical result. Because the spacing between the four levels of educational experience is very uneven, the meaning of this average would be very questionable. In short, an average requires a variable to be interval. Sometimes you have variables that are "in between" ordinal and interval, for example, a five-point likert scale with values "strongly agree", "agree", "neutral", "disagree" and "strongly disagree". If we cannot be sure that the intervals between each of these five values are the same, then we would not be able to say that this is an interval variable, but we would say that it is an ordinal variable. However, in order to be able to use statistics that assume the variable is interval, we will assume that the intervals are equally spaced.

Does it matter if my dependent variable is normally distributed?

When you are doing a t-test or ANOVA, the assumption is that the distribution of the sample means are normally distributed. One way to guarantee this is for the distribution of the individual observations from the sample to be normal. However, even if the distribution of the individual observations is not normal, the distribution of the sample means will be normally distributed if your sample size is about 30 or larger. This is due to the "central limit theorem" that shows that even when a population is non-normally distributed, the distribution of the "sample means" will be normally distributed when the sample size is 30 or more, for example see "Central limit theorem demonstration."http://www.ats.ucla.edu/STAT/stata/ado/teach/clt.htm

If you are doing a regression analysis, then the assumption is that your residuals are normally distributed. One way to make it very likely to have normal residuals is to have a dependent variable that is normally distributed and predictors that are all normally distributed, however this is not necessary for your residuals to be normally distributed. You can see Regression with SAS: Chapter 2 - Regression Diagnostics, Regression with SAS: Chapter 2 - Regression Diagnostics, or Regression with SAS: Chapter 2 - Regression Diagnostics.

The Statistic Analytic Reporter

1. Overview

The Enterprise line of STATISTICA products is designed for multi-user, collaborative analysis applications. STATISTICA analytic procedures are deployed on a server and are architected for the simultaneous use by any number (tens, hundreds, thousands) of users. An important component of the STATISTICA Enterprise architecture is the STATISTICA Analytic Reporter. The STATISTICA Analytic Reporter allows report configurations and report generation in STATISTICA Enterprise to deliver formatted documents (PDF, HTML, MS Word) and analysis summaries of any of the tabular and graphical results produced by STATISTICA.

2. Business Applications

The STATISTICA Enterprise software is used in every industry for countless applications so a comprehensive overview of the myriad ways in which the STATISTICA Analytic Reporter is used is not possible. In short, the STATISTICA Analytic Reporter is used for the following categories of applications:

2.1. Production of Standard, Formatted Reports

In a variety of industries, the documentation of ongoing analysis summaries for decision-making and historical record is a mandated aspect of doing business. For example, in FDA-regulated industries, standard reports are inherent to Good Manufacturing Practices (GMP) to demonstrate ongoing product quality monitoring and due diligence in determining potential trends and root causes. These Annual or Periodic Product Reviews (APRs/PPRs) mandated for GMP compliance can be time-consuming to produce. The STATISTICA Analytic Reporter streamlines the time to produce APRs/PPRs, allowing Quality Control personnel to focus on the review of the historical charts, and reduces the potential for errors in the analyses in the charts by automating the process.

2.2. Delivery of Snapshots and Updates to Consumers of Information

In any Enterprise, there are active, "power" users of STATISTICA, and the periodic consumers of information whose responsibilities and needs require that they are informed of the "big picture." For the consumers of information, who need to always be aware of "what happened yesterday", batch report summaries provide a convenient mechanism to stay informed. The report summaries can be waiting for the user when the user logs into the STATISTICA Knowledge Portal, posted to an Intranet site where the user can browse as needed, or can be automatically printed. For example, one of StatSoft's customers automatically prints a daily report every morning at 6AM for pickup by a key consumer for their review when they arrive at 7AM. The reports are in a standard format, so the review is streamlined to allow for the selection of focus on what the important issues of the day are.

2.3. Dashboards One output of the STATISTICA

Analytic Reporter are HTML documents with graphical results as JPEGs/PNGs/GIFs. The HTML documents outputted by the STATISTICA Analytic Reporter are often the sub-components feeding into a dashboard summary of Key Performance Indicators (KPIs) within an organization. The automatic scheduling of updates to the standard reports allow for the graphical summaries of KPIs to be refreshed at the desired frequency to suit the business need.

3. How It Works

Below is a summary of the details of the STATISTICA Analytic Reporter, and the flexibility it provides as the report templating and generation mechanism within the STATISTICA Enterprise platform.

3.1. Analysis Automation The STATISTICA

Enterprise application platform includes the industry-standard STATISTICA Visual Basic (SVB) for analysis automation. The complete set of STATISTICA analyses and graphics are exposed for automation and macro recording. This analysis automation using STATISTICA Visual Basic is the first step in configuring a STATISTICA Report configuration.

3.2. Report Templating

The STATISTICA Analytic Reporter provides a graphical user interface for defining the layout of formatted documents, including the placement of graphs and tables, the contents and formatting of headers/footer, static and dynamic text elements, and any additional formatting elements specific to the document type. The results of the report template definition are saved as a STATISTICA Report template document. Importantly, these Report Templates are stored centrally, wrapped with user access control and security, in the STATISTICA Configurations Database, deployed on an industry-standard relational database management system (RDBMS).

3.3. Report Generation Engine

At the heart of the STATISTICA Analytic Reporter is the report generation engine. Its role is to retrieve the latest version of the respective report template, "pickup" the latest tabular and graphical results for the associated analyses, populate the report template with the latest results, and then lastly generate the output document in the desired format. The STATISTICA Analytic Reporter handles all of these steps automatically.

3.4. Scheduling and Batching

The STATISTICA Analytic Reporter may either be run on-demand or as batch, scheduled tasks. Both approaches have their place for respective applications within an Enterprise. With on-demand report, users may be prompted to input application options. For example, with Annual Product Review reports, the user may be prompted to select a Product and a Date Range. For batch applications, the STATISTICA Analytic Reporter has its own batch scheduler for configuring the timing and frequency for each individual report. The STATISTICA Analytic Reporter also supports batch scheduling via the Windows AT Scheduler.

3.5. Reports Delivery

The STATISTICA Analytic Reporter delivers updated reports documents to a designated. By default, the location is the STATISTICA Knowledge Portal, a secure document repository integrated into STATISTICA Enterprise's access control and security infrastructure. The STATISTICA Knowledge Portal allows for the organized, browsing of the latest reports to which a respective user has permissions.