Graphics:

### Normal distribution
curve(dnorm(x),xlim=c(-4,4),ylab="density",main="Normal Distribution")

# 2 normal curves
curve(dnorm(x,mean=0,sd=5), xlim=c(-20,50), col="red")
curve(dnorm(x,mean=25,sd=5), add=TRUE, col="blue")


### y^2 = 4*x
y = seq(0, 5, by=0.01)
y.upper = 2*sqrt(x)
y.lower = -2*sqrt(x)
y.max = max(y.upper)
y.min = min(y.lower)
plot(c(-2,5), c(y.min,y.max), type="n", xlab="x", ylab="y")
lines(x, y.upper)
lines(x, y.lower)
abline(v=-1)
points(1,0)
text(1,0,"focus (1,0)", pos=4)
text(-1,y.min,"directrix x = -1", pos=4)
title("The parabola y^2 = 4*x")

### x^2 - y^2/3 = 1
plot(c(-max(x),max(x)),c(-max(y),max(y)),type='n',xlab='x',ylab='y')
x = seq(1,5,by=0.01)
y = sqrt(3*x^2-3)
lines(x,y); lines(-x,y); lines(x,-y); lines(-x,-y);
abline(0,sqrt(3)); abline(0,-sqrt(3))
points(2,0); points(-2,0)
text(2,0,"focus (2,0)",pos=4)
text(-2,0,"focus (2,0)",pos=2)
text(5,8.5,"asymptote y = sqrt(3)*x", pos=2)

### 
opar1 = par(las=1, mar=c(4,4,3,2))
plot(cars$speed,cars$dist,axes=FALSE,xlab="",ylab="",type="n")
points(cars$speed, cars$dist,
       col=ifelse(cars$speed>9,"darkseagreen4","red"),
       pch=ifelse(cars$speed>9,1,3))
axis(1);axis(2)

opars = par(las=0)
mtext("Speed (mph)",side=1,line=3)
mtext("Dist (ft)",side=2,line=3)

box()

legend(x=10,y=100,c("Low","High"),
       col=c("darkseagreen3","red"), 
       pch=c(3,1),bty="n")

### 
curve(100*(x^3-x^2)+15, from=0, to=1,
      xlab=expression(alpha),
      ylab=expression(100 %*% (alpha^3 - alpha^2) +15),
      main=expression(paste("Function : ", 
           f(alpha) == 100 %*% (alpha^3 - alpha^2) + 15)))
myMu = 0.5; mySigma = 0.25
par(usr=c(0,1,0,1))
text(0.1,0.1,bquote(sigma[alpha] == .(mySigma)), cex=1.25)
text(0.6,0.6,paste("(The mean is ", myMu, ")", sep=""), cex=1.25)
text(0.5,0.9,bquote(paste("sigma^2 = ", sigma^2 == .(format(mySigma^2,2)))))

GEE

### GEE ?
library(nlme)
lme(travel~1, random=~1|Rail, method="ML", data=Rail)

y<-Rail$travel
id<-Rail$Rail
time<-rep(1:3,6)

X<-rep(1,length(y))
mu<-solve(t(X)%*%X)%*%t(X)%*%y

Y<-cbind(y[time==1],y[time==2],y[time==3])
Res<-Y-matrix(rep(mu,18),6)
R<-cor(Res)
phi<-1/var(y)

V<-matrix(rep(1,9),3)
for(i in 1:3) for(j in 1:3) if(i!=j) V[i,j]=mean(c(R[1,2],R[1,3],R[2,3]))

solve(t(X)%*%solve(kronecker(V,diag(6)))%*%X)%*%t(X)%*%solve(kronecker(V,diag(6)))%*%y

alpha<-mean((Res[,1]*Res[,2]+Res[,2]*Res[,3])/2)/phi

Logistic Regression - estimation

### GLM

## logistic regression
respire<-read.csv("F:/publish/data analysis using R/data/respire.csv")
respire2<-read.csv("F:/publish/data analysis using R/data/respire2.csv")

out<-glm(outcome ~ treat, weights=count, family=binomial, data=respire)

logist.lm<-function(beta){
p0<-exp(beta[1])/(1+exp(beta[1]))
p1<-exp(beta[1]+beta[2])/(1+exp(beta[1]+beta[2]))
ll<-16*log(p0)+48*log(1-p0)+40*log(p1)+20*log(1-p1)
return(-ll)
}
optim(c(-1,2),logist.lm)

library(rootSolve)

model<-function(beta){
x<-(respire$treat=="test")
y<-respire$outcome
mu<-exp(beta[1]+beta[2]*x)/(1+exp(beta[1]+beta[2]*x))
W<-respire$count
F1<-sum(W*(y-mu))
F2<-sum(W*(y-mu)*x)
c(F1=F1,F2=F2)
}
multiroot(f=model,start=c(-1,2))$root

Mixed Models - REML

library(nlme)

# 6 Rails have 3 repeatitions each.
n<-6; r<-3;

# -log Restricted Maximum Likelihood
# beta[1]: mean, exp(beta[2]): sigma_a, exp(beta[3]): sigma
# exp is used to prevent negative estimator for variance
reml<-function(beta){
y<-Rail$travel
X<-rep(1,n*r)
Z<-kronecker(diag(n),rep(1,r))
V<-exp(beta[2])*Z%*%t(Z)+exp(beta[3])*diag(n*r)

ml<-log(det(V))/2+t(y-beta[1])%*%solve(V)%*%(y-beta[1])/2+log(det(t(X)%*%solve(V)%*%X))/2

return(ml)
}

# minimize reml function with optim()
optim(c(66.5, 5, 2),reml)

$par
[1] 66.505645 6.421954 2.782963
> sqrt(exp(c(6.421954,2.782963)))
[1] 24.803307 4.020802

# ML: 22.62435 4.020779
# REML: 24.80547 4.020779
lme(travel~1, random=~1|Rail, data=Rail, method="ML")
lme(travel~1, random=~1|Rail, data=Rail)

Simulation

### Functions
test_stat <- function(x){
  ts <- numeric()
  
  X <- matrix(x,r)
  rs <- rowSums(X)  # ni.
  cs <- colSums(X)  # n.j
  
  # E(n)  
  M <- outer(rs,cs)/sum(X)  
    
  # 1st test statistic (LRT)
  ts[1]<-1-pchisq(2*sum(X*log(X/M),na.rm=T),(r-1)*(c-1));

  # 2nd test statistic (Pearson's Chi-square)
  ts[2]<- 1-pchisq(sum((X-M)^2/M),(r-1)*(c-1));
  
  ts[3:12]<-exact(X,rs,cs);
  
  return(ts)
}
exact <- function(x,rs,cs){
  T <- r2dtable(iter,rs,cs)
  obs <- unlist(mmh(x))
  exp <- sapply(T,mmh)

  return(apply(exp>obs,1,mean))
}
mh<-function(x,i,j){
  N<-x[1,1]+x[1,j]+x[i,1]+x[i,j];
  return((x[i,j]*x[1,1]/N-x[i,1]*x[1,j]/N)/sqrt((x[1,1]+x[i,1])/N*(x[1,1]+x[1,j])/N*(x[i,j]+x[i,1])/(N-1)*(x[i,j]+x[1,j])))
}
mmh<-function(x){
  Z1<-mh(x,2,2);  Z2<-mh(x,2,3);  Z3<-mh(x,3,2);  Z4<-mh(x,3,3)
  Z_one<-max(Z1,Z2,Z3,Z4)
  Z_two<-max(abs(Z1),abs(Z2),abs(Z3),abs(Z4))
  return(list(Z_one,Z1,Z2,Z3,Z4,
              Z_two,abs(Z1),abs(Z2),abs(Z3),abs(Z4)))
}


### sample size and number of iteration
n <- 100; r<-3; c<-3; s <- 1000; iter <- 1000; 

### test statistics
num_ts<-12; names_ts <- c('QL', 'Qp',
                          'One-sided', 'MH22', 'MH23', 'MH32','MH33',
                          'Two-sided', 'MH22', 'MH23', 'MH32','MH33')

### parameters
alpha <- c(1,0.2,0.2); beta <- c(1,0.2,0.1); gamma <- matrix(c(1,1,1,1,3.,1.,1,1.,1.),r)

### s 2*2 random tables with n from multinomial distribution
tables <- rmultinom(s, n, outer(alpha,beta)*gamma)

p <- apply(tables,2,test_stat)
rownames(p)<-names_ts

### empirical tail probability
percent <- cbind(apply((p<=0.05),1,mean,na.rm=TRUE))
rownames(percent)<- names_ts

write.csv(rbind(alpha,beta,gamma),"E:/documents/maxor/out/sim.csv",append=TRUE)
write.csv(percent,"E:/documents/maxor/out/sim.csv",append=TRUE)

Nonlnear Models - GLS

### Indometh

ols<-nls(conc~SSbiexp(time,b1,b2,b3,b4),
              data=Indometh[Indometh$Subject==1,])
beta<-coef(ols)

for(i in 1:10){
  oldbeta<-beta
  out<-nls(conc~SSbiexp(time,b1,b2,b3,b4),
             data=Indometh[Indometh$Subject==1,],
             weights=1/SSbiexp(time,beta[1],beta[2],beta[3],beta[4])^2)
  beta<-coef(out)              
}

BETA<-matrix(rep(0,6*4),6)
for(i in 1:6){
  BETA[i,]<-coef(nls(conc~SSbiexp(time,b1,b2,b3,b4),
            data=Indometh[Indometh$Subject==i,]))
}
colnames(BETA)<-c('b1','b2','b3','b4')


### Theoph
coplot(conc~Time | Subject, data=Theoph)
nls(conc~SSfol(Dose,Time,lKe,lKa,lCl),data=Theoph)
ols<-nls(conc~SSfol(Dose,Time,lKe,lKa,lCl),
              data=Theoph[Theoph$Subject==1,])
beta<-coef(ols)

BETA<-matrix(rep(0,12*3),12)
for(i in 1:max(as.numeric(Theoph$Subject)))
  BETA[i,]<-coef(nls(conc~SSfol(Dose,Time,b1,b2,b3),
                 data=Theoph[Theoph$Subject==i,]))

Poisson Regression

melanoma<-read.csv("U:/data/R/datasets/melanoma.csv")

out<-glm(cases~age+region, family=poisson(link="log"), 
                           offset=log(total), data=melanoma)
out1<-glm(cases~age, family=poisson(link="log"), 
                     offset=log(total), data=melanoma)
out2<-glm(cases~region, family=poisson(link="log"), 
                        offset=log(total), data=melanoma)

# type 3
anova(out1,out)
anova(out2,out)

# melanoma
Obs     age     region    cases     total  

  1    35-44    south       75      220407 
  2    45-54    south       68      198119 
  3    55-64    south       63      134084 
  4    65-74    south       45       70708 
  5    75+      south       27       34233 
  6    <35      south       64     1074246 
  7    35-44    north       76      564535 
  8    45-54    north       98      592983 
  9    55-64    north      104      450740 
 10    65-74    north       63      270908 
 11    75+      north       80      161850 
 12    <35      north       61     2880262 

Datasets

# PK
xyplot(conc~time|Subject, data=Indometh) # IV
xyplot(conc~Time|Subject, data=Theoph) # oral

# Linear
plot(dist~speed, data=cars)
plot(eruptions~waiting, data=faithful)
xyplot(weight~Time | Diet, data=ChickWeight)
lme(weight~Time+Diet, random=~1+Time |Chick, data=ChickWeight)
xyplot(height~age|Seed, data=Loblolly)

# Latin-square
OrchardSprays

# ANOVA
boxplot(weight~group, data=PlantGrowth)
boxplot(weight~feed, data=chickwts)

boxplot(breaks~wool+tension, data=warpbreaks)


# Nonlinear
xyplot(uptake~conc|Plant, data=CO2)
xyplot(density ~ conc | Run, data=DNase)
xyplot(rate~conc|state, data=Puromycin)
xyplot(circumference~age|Tree, data=Orange)
plot(pressure~temperature, data=pressure)

# multivariate
attitude
trees

# logistic regression
esoph
glm(cbind(ncases, ncontrols)~agegp+tobgp*alcgp,data=esoph,family=binomial())

# paired t-test
t.test(extra~group, paired=TRUE, data=sleep)

##### MASS
# Survival
Melanoma
VA
gehan
leuk

xyplot(BPchange~Dose | Animal+Treatment, data=Rabbit)
xyplot(size~Time | treat, data=Sitka)

# logistic regression
birthwt

# Poisson Regression
boxplot(log(y)~limit, data=Traffic)
summary(glm(y~limit, data=Traffic, family=poisson))
summary(glm(y~day+limit, data=Traffic, family=poisson))

# nonlinear
plot(steam)

# two-way anova
genotype

# ANACOVA
anova(lm(Postwt~Prewt+Treat, data=anorexia))

# paired t-test
write.csv(anorexia[anorexia$Treat=="CBT"2:3],"u:/data/R/datasets/paired_CBT.csv",row.names=FALSE)

McNemar's test

# simulate a dataset
group<-rep(c('case','cont'),rep(100,2))
pre<-rbinom(200,1,0.7)
post<-rbinom(200,1,0.8)
mcnemar.test(pre,post)

# create a variable post - pre
matched<-data.frame(pre,post,group)
matched$change<-post-pre

# McNemar's test from logistic regression (0 has to be deleted)
summary(glm(factor(change)~1, data=matched[matched$change!=0,],family=binomial))

summary(glm(factor(change)~group, data=matched[matched$change!=0,],family=binomial))

Multiple Comparisons

### Use warpbreaks dataset
amod<-aov(breaks~tension, data=warpbreaks)

### 1-way ANOVA (exact)
library(multcomp)
summary(glht(amod, linfct=mcp(tension="Dunnett")))
summary(glht(amod, linfct=mcp(tension="Tukey")))

# CI
confint(glht(amod, linfct=mcp(tension="Dunnett")), level=0.9)

# create a contrast
contr<-rbind("M-L"=c(-1,1,0), "H-L"=c(-1,0,1), "H-M"=c(0,-1,1))
glht(amod, linfct=mcp(tension=contr))

# Tukey test
TukeyHSD(chickwts.out,"feed")

# Closed test
pairwise.t.test(warpbreaks$breaks, warpbreaks$tension, p.adj="none")
pairwise.t.test(warpbreaks$breaks, warpbreaks$tension, p.adj="bonf")
pairwise.t.test(warpbreaks$breaks, warpbreaks$tension, p.adj="hommel")

ROC curve

library(ROC)

R1 <- rocdemo.sca( rbinom(40,1,.3), rnorm(40), caseLabel="new case", markerLabel="demo Marker" );

plot(R1, line=TRUE, show.thresh=TRUE);

plot(R1, line=TRUE);

reshape

# wide to long
d1<-data.frame(subject=c("an1","an2"),eat1=1:2,eat3=5:6,trt=c("t1","t2"))
d1.1<-reshape(d1, idvar="subject",varying=list(c("eat1","eat3")),v.names="eat",direction="long")
d1.2<-reshape(d1, idvar="subject",varying=list(c("eat1","eat3")),v.names="eat",times=c(1,3),timevar="week",direction="long")

d2<-data.frame(subject=c("an1","an2"),eat1=1:2,eat3=5:6,gain1=14:15,gain3=20:21,trt=c("t1","t2"))
d2.1<-reshape(d2,idvar="subject",varying=list(c("eat1","eat3"),c("gain1","gain3")),v.names=c("eat","gain"),times=c(1,3),timevar="week",direction="long")
d2.2<-reshape(d2,idvar="subject",varying=list(c("eat1","eat3")),v.names="eat",times=c(1,3),timevar="week",direction="long")
d2.3<-reshape(d2,idvar="subject",varying=list(c("eat1","eat3")),drop=c("gain1","gain3"),v.names="eat",times=c(1,3),timevar="week",direction="long")

# long to wide
reshape(d1.2,idvar="subject",timevar="week",v.names="eat",direction="wide")
reshape(d2.1,idvar="subject",timevar="week",v.names=c("eat","gain"),direction="wide")

from Reshaping data with the reshape function in R by Soren Hojsgaard

Paired t-test with Mixed Models

library(MASS)
library(nlme)

# simulated a dataset
paired<-mvrnorm(n=30,mu=c(10,10.2),Sigma=matrix(c(5,4,4,5),2,2))

before<-paired[,1]
after<-paired[,2]


# paired t-test
t.test(before,after,paired=TRUE)



### Mixed Models for a paired t-test 
# reshape data for mixed models
y<-c(before,after)
id<-rep(c(1:length(before)),2)
time<-c(rep(0,length(before)),rep(1,length(after)))


# fit mixed models
summary(lme(y~time, random=~1|id))

Mixed Models - longitudinal data

### Sitka (79 trees with 5 time points)

# random intercept models = compound symmetry
sitka.lme <- lme(size~treat*ordered(Time), random=~1|tree, data=Sitka)

Mixed Models - random coefficient models

library(nlme)

coplot(Y ~ EP|No, data=petrol)

Petrol <- petrol
Petrol[,2:5] <- scale(Petrol[,2:5],scale=F)

# Y~EP for each No; no overal mean
pet1.lm <- lm(Y~No/EP-1, data=Petrol)

# Y~ mu_i + beta*EP; no overal mean
pet2.lm <- lm(Y~No-1+EP, data=Petrol)
pet3.lm <- lm(Y~SG+VP+V10+EP, data=Petrol)

# Random intercept model
pet3.lme <- lme(Y~SG+VP+V10+EP, random=~1|No, data=Petrol)

# To compare likelihood, use "ML"
pet3.lme <- update(pet3.lme, method="ML")
anova(pet3.lme, pet3.lm)
pet4.lme <- update(pet3.lme, fixed=Y~V10+EP)
# pet4.lme <- lme(Y~V10+EP, random=~1|No, method="ML", data=Petrol)
anova(pet3.lme, pet4.lme)

## Random intercept + Random coefficient
pet5.lme <- lme(Y~V10+EP, random=~1+EP|No, method="ML", data=Petrol)
anova(pet4.lme, pet5.lme)

Mixed Models - RCBD

### Randomized Complete Block Design
library(nlme)

rcb <- read.csv("rcb.csv")
rcb$ingot <- as.factor(rcb$ingot)
rcb$metal <- as.factor(rcb$metal)

### lme
rcbFit <- lme(pres~ metal, random=~1|ingot, data=rcb)

### manually
t<-3; b<-7

mean_b<-tapply(rcb[,3],rcb[,1],mean)
mean_t<-tapply(rcb[,3],rcb[,2],mean)[t:1]

msa<-b*sum((mean_t-mean(mean_t))^2)/(t-1)
msb<-a*sum((mean_b-mean(mean_b))^2)/(b-1)
mse<-sum((y-rep(mean_t,b)-rep(mean_b,rep(t,b))+mean(y))^2)/(t-1)/(b-1)
F<-msa/mse
p_value <- 1-pf(F,(t-1),(t-1)*(b-1))

sb<-(msb-mse)/t # sigma_b

/* data & SAS code */
data rcb;
   input ingot metal $ pres @@;
   datalines;
1 n 67.0  1 i 71.9  1 c 72.2
2 n 67.5  2 i 68.8  2 c 66.4
3 n 76.0  3 i 82.6  3 c 74.5
4 n 72.7  4 i 78.1  4 c 67.3
5 n 73.1  5 i 74.2  5 c 73.2
6 n 65.8  6 i 70.8  6 c 68.7
7 n 75.6  7 i 84.9  7 c 69.0
run;

proc mixed data=rcb;
   class ingot metal;
   model pres=metal;
   random ingot;
run;

Structural Equation Model

### Principle Component Analysis

## with prcomp()
prcomp(USArrests, scale=T)
prcomp(~Murder+Assault+Rape, data=USArrests, scale=T)
plot(prcomp(~Murder+Assault+Rape, data=USArrests, scale=T))
summary(prcomp(~Murder+Assault+Rape, data=USArrests, scale=T))
biplot(prcomp(~Murder+Assault+Rape, data=USArrests, scale=T))

## with princomp()
princomp(USArrests, cor=T)
princomp(~Murder+Assault+Rape, data=USArrests, cor=T)
summary(princomp(~Murder+Assault+Rape, data=USArrests, cor=T))
princomp(~Murder+Assault+Rape, data=USArrests, cor=T)$score
plot(princomp(~Murder+Assault+Rape, data=USArrests, cor=T))
biplot(princomp(~Murder+Assault+Rape, data=USArrests, cor=T))

# missing handling
USArrests[1,2]<-NA
princomp(~Murder+Assault+Rape, data=USArrests, na.action=na.exclude, cor=T)


### Factor Analysis with factanal()
factanal(~Murder+Assault+UrbanPop+Rape, factors=1, data=USArrests)

### FactoMineR

Collapse table and calculate p-value by Monte Carol Method

# colapse rows and columns whose cell counts < 4
rcbind <- function(tt)
{
  # rearraign table
  cc <- names(sort(tt[names(sort(rowSums(tt==max(tt)), decreasing=T))[1],], decreasing=T))
  rr <- names(sort(tt[,names(sort(colSums(tt==max(tt)), decreasing=T))[1]], decreasing=T))
  tt <- tt[rr,cc]

  # p-value for Consensus Test (binomial test)
  m1=max(tt)
  m2=max(tt[1,2:dim(tt)[2]],tt[2:dim(tt)[1],1])
  p <- binom.test(m1,m1+m2)$p.value # 2-sided

  if(p>=0.05){
    print("There is no concensus.")
    return(list(0,tt))
  } else
  if(sum(tt[2:dim(tt)[1],2:dim(tt)[2]]>=4)==0){
    print("There is no double mutation cell >=4.")
    return(list(0,tt))
  } else {
    # inner table
    inner <- tt[2:dim(tt)[1],2:dim(tt)[2]]
    rrr <-rowSums(inner>=4)
    ccc <-colSums(inner>=4)

    # colnames and rownames whose cell counts are larger than 4
    rnames <- c(rownames(tt)[1],names(rrr[rrr>=1]))
    cnames <- c(colnames(tt)[1],names(ccc[ccc>=1]))

    # colnames and rownames which is collapsed
    mr <- rownames(tt)[-charmatch(rnames,rownames(tt))]
    mc <- colnames(tt)[-charmatch(cnames,colnames(tt))]

    # colapsing rows and columns whose cell counts are smaller than 4
    return(list(1,rbind(cbind(tt[rnames,cnames],rowSums(tt[rnames,mc]))
          ,c(colSums(tt[mr,cnames]),sum(tt[mr,mc])))))
  }
}

test_stat <- function(x){
  # Pearson's Chi-square test
  expect <- (matrix(rowSums(x))%*%t(matrix(colSums(x))))/sum(x)
  x2 <- sum((x-expect)^2/expect)

  # LRT
  temp <- x*log(x/expect)
  g2 <- 2*sum(temp[x!=0])

  # Maximal Odds Ratio & Maximal MH
  or <- rep(0,(dim(x)[1]-1)*(dim(x)[2]-1))
  dim(or) <- c(dim(x)[1]-1,dim(x)[2]-1)
  diff <- or
  for(i in 2:dim(x)[1])
    for(j in 2:dim(x)[2]){
       or[i-1,j-1] <- abs(log(x[i,j]*x[1,1]/x[i,1]/x[1,j])
                    /sqrt(1/x[1,1]+1/x[1,j]+1/x[i,1]+1/x[i,j]))
       diff[i-1,j-1] <- abs(x[i,j]*x[1,1]-x[i,1]*x[1,j])
                 /sqrt((x[1,1]+x[i,1])*(x[1,1]+x[1,j])*(x[1,j]+x[i,j])
                 *(x[i,1]+x[i,j])/(x[i,j]+x[1,1]+x[i,1]+x[1,j]-1))

    }
#  df <- (dim(x)[1]-1)*(dim(x)[2]-1)
#  1-pchisq(c(x2,g2),df) 
  
  c(x2,g2,max(or),max(diff))

}

 

exact <- function(x){

  ts_o <- test_stat(x)

  reference <- r2dtable(iter, rowSums(x), colSums(x))

  ts_e <- sapply(reference, test_stat)

  p_value <- c(fisher.test(x)$p.value,apply(ts_e>ts_o,1,mean,na.rm=T))
  names(p_value) <- c("Exact","X2","LRT","Maximal OR","Maximal MH")
  
  return(p_value)

}

nsi<-read.csv("g:/programming/R/nsi.csv")

iter =10000
position1 <- nsi$x18
position2 <- nsi$x19

ss <- rcbind(table(position1,position2))

if(ss[[1]]) p_value <- exact(ss[[2]])

ss[[2]]
p_value

aa <- r2dtable(1,c(90,50,40),c(90,45,45))
exact(aa[[1]])

Mortgage

mortgage <- function(P, r, year, py){
  m <- year*12;
  Ei <- rep(0,m);
  mr <- r/12;
  E <- P*mr/((1+mr)^m -1);
  pmt <- P*mr+E;

  for(i in 1:m) Ei[i]<-E*(1+mr)^(i-1);
  Ii <- pmt-Ei;

  cat(year, "year at", r*100, "%\n");
  cat("Monthly Payment:", pmt, "\n");
  cat("Total Payment:", pmt*m, ", Total Interest:", pmt*m-P, "\n");
  cat("Total Principal:", sum(Ei[1:(12*py)]), 
      ", Total Interest:", sum(Ii[1:(12*py)]),"(After", py, "years)\n");
}

P <- 210000 # principal
r <- 0.05625 # interest rate
year <- 15 # loan period
py <- 5 # short term

mortgage(P,r,year,py)