Simpson paradox_Longitudinal study

x1 <- c(1,2,3,4)
y1 <- x1 + 5

x2 <- x1 + 7
y2 <- x2 – 7

x <- c(x1,x2)
y <- c(y1,y2)

postscript(“simpson.eps”, paper=”special”, width=4.5, height=3)

par(las=1)
par(mar=c(3,3,0.5,0.5))
par(mgp=c(2,1,0))

plot(x,y, cex=2, pch=21,
col=rep(c(“blue”, “red”), each=4), bg=rep(c(“lightblue”, “pink”), each=4),
xlim=range(x)+c(-2,2), ylim=range(y)+c(-2,2))
abline(lm(y1 ~ x1), col=”blue”, lwd=2)
abline(lm(y2 ~ x2), col=”red”, lwd=2)
abline(lm(y ~ x), lwd=2, lty=2)

dev.off()

freq table

require(gmodels)data(infert, package = "datasets")CrossTable(infert$education, infert$induced, expected = TRUE)CrossTable(infert$education, infert$induced, expected = TRUE, format="SAS")CrossTable(infert$education, infert$induced, expected = TRUE, format="SPSS")CrossTable(warpbreaks$wool, warpbreaks$tension, dnn = c("Wool", "Tension"))

function findlink(pkg, fn) { var Y, link; Y = location.href.lastIndexOf(“\\”) + 1; link = location.href.substring(0, Y); link = link + “../../” + pkg + “/chtml/” + pkg + “.chm::/” + fn; location.href = link; }

increas the # of tick intervals

library(plotrix)
hist (data1, seq(-0.01,0.3,0.01),freq = FALSE, main=”Figure1.”,xlab=”return”,lab=c(1,1,12),xaxt=”n”) staxlab(1,at=seq(-0.01,0.3,by=0.02),labels=seq(-0.01,0.3,by=0.02))

barplot(one,axes=F)
axis(2,at=seq(0,80,by=5))

descriptive statistics

require(psych)
describe(x)
describe.by(x,group)

It gaves
A data.frame of the relevant statistics:
item name
item number
number of valid cases
mean
standard deviation
median
mad: median absolute deviation (from the median)
minimum
maximum
skew
kurtosis
standard error

statmod <- function(x) {
z <- table(as.vector(x))
names(z)[z == max(z)]
}

Demo R

\documentclass[hyperref={dvips,colorlink,bookmarks}]{beamer}

\mode {
\definecolor{one}{rgb}{0.2,0.5,0.9} %0.1
\usecolortheme[named=one]{structure}
\usetheme{Madrid} %CambridgeUS Szeged
\useinnertheme{rectangles}
%\usecolortheme{rose}
\usepackage{multimedia}
\usepackage{amssymb,graphicx,listings}
\useoutertheme{infolines} %miniframes infolines
}

\title{Demo R}
\author{Advanced stat TA}
\institute{NTU AGEC}
\date{\today}

\begin{document}
\lstset{numbers=left,tabsize=5,numberstyle=\tiny,
keywordstyle=\color{blue!80},
commentstyle=\color{black!40!green!100},
basicstyle=\scriptsize,stepnumber=1,numbersep=5pt, language={R}}

\frame{\titlepage}

\frame{\frametitle{Demo 1}
\begin{block}{Ex1.1-4}
EX1.1-3 The $National \; geographic$
reported the following litter sizes for 20 lions.
\end{block}

\begin{itemize}
\item Construct a freq table. (\textcolor{blue}{stem() ,table(),length()})
\item Draw a freq histogram. (\textcolor{blue}{barplot()} )
\item Draw a relative freq histogram.
\item Is there a typical litter size.
\end{itemize}
Making yourself R

\begin{enumerate} []
\item locate your data file. (\textcolor{blue}{setwd()})
\item import data from the working space.(\textcolor{blue}{scan()})
\item coding…
\end{enumerate}
}

\begin{frame}[fragile]
\frametitle{Demo1 code}
\begin{block}{}
\begin{lstlisting}

setwd("c:/data/advsta")
one<-scan("Exercise 1_1-03.txt",sep=",")

#Q1-1
stem(one)
t1=table(one)
t2=t1/length(one)

#Q1-2~4
par(mfrow=c(1,2),bg="white",col.axis="black")
barplot(t1,ylab="Freq")
text(2,9,"<=typical little size",col="blue",cex=1,adj = c(-0.4,0))
barplot(t2,ylab="Density",col="lightblue")

\end{lstlisting}
\end{block}
\end{frame}

\begin{frame}
\frametitle{Demo2}
\begin{block}{Ex1.1-4}
Some ornithologists were interested in the clutch size of the common gallinule.
Thye observed the number of eggs in each of 117 nests.
\end{block}
\begin{itemize}
\item Construct a freq table.
\item Draw a histogram.
\item What is the mode?
\end{itemize}

\end{frame}
\end{document}

Three Asymptotic Tests: LR, Wald, LM test

LR test : restricted parameter vector and specify distribution(∵ likelihood function)
Wald test : unrestricted parameter vector only
LM test: restricted parameter vector only, the popular method in econometrics since restricted model is easy to get.

Method 1

#read data
hs1 <- read.table("http://www.ats.ucla.edu/stat/R/notes/hs1.csv", header=T, sep=",") attach(hs1) # classify variable honcomp = 60
honcomp[1:20]
#LR test for logistic reg
lr <- glm(honcomp~female+read, family=binomial) lr2 <- glm(honcomp~female+read+science+math+as.factor(prog), binomial, na.action=na.exclude) drop1(lr2, test="Chisq")

Method 2

see http://cran.r-project.org/src/contrib/Descriptions/lmtest.html
#download the package lmtest
library(lmtest)
lrtest(lr,lr2)
waldtest(lr,lr2)
lmtest(lr,lr2)


Optim method

大森裕浩 本文

見本用に関数最適化ルーチン “optim” の html ヘルプファイルを別途置いておきます。 “optim” はオプションで複数の最適化手法を選ぶことができ、問題に応じた適切な手法を簡単に試行錯誤で探すことができます。また複数の手法で得られた解を比較すれば、解の妥当性を容易に検証することもできます。これだけでも R を使う価値があると思います。 いずれも複数パラメータが可で、グラディエントベクトルを明示的にあたえたり、数値微分で計算するように指示できます。otptim は再帰的に使用できます(私の経験では、制御パラメータ値の選択に極めて敏感な、たちの悪い最適化問題で、制御パラメータ値自身を変数に取る二段階の最適 化で良い結果を得ることが出来ました)。

===================================================================
(1) Nelder-Mead法 関数値のみを使用、頑健、相対的に遅い(既定手法) 。
———————————————————————————————–
(2) BFGS法(準ニュートン法) 関数値とグラディエントを使用 。
———————————————————————————————–
(3) CG法(共役勾配法) BFGS法より破綻しやすいが、メモリ使用量が少く大規模問題に
適する三種類の更新法(Fletcher-Reeves, Polak-Ribiere, Beale-Sorenson)を選択可 。
———————————————————————————————–
(4) L-BFGS-B法 矩形型拘束条件を持つ準ニュートン法 。
———————————————————————————————
(5) SANN法 Metropolis 法を使うシミュレーテッドアニーリング法、関数値のみを使用、
頑健、遅い 。
==================================================================

9*9 table in R and c++

How to use R write 9*9 table

R code

for(i in seq(2,9,by=1)){
for(j in 1:9){
cat(i,"*",j,"=",i*j,"\n")
}}


Cpp

#include
using namespace std;
int main() {
for(int j = 1; j = 10)
cout << "=" <
using namespace std;
int main() {
for (int i = 2, j = 1;
j =10)? "=": "= ") ;
if(i == 9)
cout <

Hypothesis test and distribution

A visual way to inspect hypothesis(one hypothesis one distribution)
R code

x<-seq(-10,10,length=400)
y1<-dnorm(x)
y2<-dnorm(x,m=3)
par(mar=c(5,4,2,1))
plot(x, y2, xlim=c(-3,8), type="n", xlab=quote(Z==frac(mu[1]-mu[2],
sigma/sqrt(n))), ylab="Density")
polygon(c(1.96,1.96,x[240:400],10), c(0,dnorm(1.96,m=3),y2[240:400],0),
col="lightgray", lty=0) #skyblue ,wheat ,tomato
lines(x, y2,lty=2)
abline(v=0,col="red")
lines(x, y1)
abline(v=3,lty=2,col="red")

polygon(c(-1.96,-1.96,x[161:1],-10), c(0,dnorm(-1.96,m=0), y1[161:1],0),
col="RoyalBlue", lty=0)
polygon(c(1.96, 1.96, x[240:400], 10), c(0,dnorm(1.96,m=0),
y1[240:400],0), col="RoyalBlue")
legend(4.2, .4, fill=c("lightgray","RoyalBlue"),
legend=expression(P(abs(Z)>1.96, H[1])==0.85,
P(abs(Z)>1.96,H[0])==0.05), bty="n")
text(0, .2, quote(H[0]:~~mu[1]==mu[2]))
text(3, .2, quote(H[1]:~~mu[1]==mu[2]+delta))

Powerdot slides

Powerdot,likes beamer, is a screen presentation tool for latex. I prefer Powerdot to Beamer,it is simple,and has discernible math texts and good support of pstricks. However, the Beamer text is certainly more sharp and “impressive” than Powerdot.So, the choice of them is trade-off.
Each of them has its build-in styles,but the style colors can be changed with writing some syntax,makes it more flexible ones. Here is the code for Powerdot:


\documentclass[paper=screen,
mode=present,clock,style=aggie,
display=slidesnotes]{powerdot}

\pddefinepalettes{my}{
\definecolor{pdcolor1}{rgb}{0,0,0}
\definecolor{pdcolor2}{HTML}{0033FF}
\definecolor{pdcolor3}{rgb}{.96,.94,.89}
\definecolor{pdcolor4}{rgb}{.89,.85,.69}
}

\pdsetup{
lf=Applied Microeconometrics,
rf=chang,
palette=my
}

</code
>
the bold texts are main body to change the colors of style aggie.
pdcolor1-text color
pdcolor2-section text or line color
pdcolor3-background color
pdcolor4-bottom area colors
after set-up, we call this set-up by adding the syntax "palette=my" in pdsetup.

Of course, there are many latex packages support screen presentation.you may lik to this site for detail