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
--------------------------------------------------