"backup"=backup.byday.2,
"monitor"=monitor.byday.2,
"adapt"=adapt.byday.2,
"smm"=smm.byday.2,
.id="measure")
ggplot(merged.byday.2, aes(x=Day,y=mean,group=measure,color=measure)) +
geom_point() + geom_line() + facet_wrap(~Rater,ncol=5) + theme_bw() +
xlim(0,120) + ylim(0,10) + ggtitle(label="Ratings by each team member of other members:  SIRIUS21")
ggplot(merged.byday.2, aes(x=Day,y=sd,group=measure,color=measure)) +
geom_point() + geom_line() + facet_wrap(~Rater,ncol=5) + theme_bw() +
xlim(0,120) + ylim(0,10) + ggtitle(label="Variability of ratings provided by each rater (s.d. across roles rated; SIRIUS21)")
merged.2 <- merged.byday.2 %>% group_by(measure,Rater) %>% summarize(mean=mean(mean,na.rm=T),
sd=mean(sd,na.rm=T))
labs <- c("Trust", "Team\n orientation", "Communications","Leadership","Backup\n behavior","Performance\nmonitoring","Adaptability","Shared\n mental model")
labs <- factor(labs,levels=labs)
merged.2$Measure<- labs[factor(merged.3$measure,levels=c("trust","torient","comm","ls","backup",    "monitor","adapt","smm" ))]
ggplot(merged.2,aes(x=Measure,y=mean,group=Rater,color=Rater)) +geom_line(size=1.5) +geom_point(size=3.5) + theme_minimal() +  scale_y_continuous(limits=c(0,10))+
theme(axis.text.x=element_text(size=12)) +
labs(title="Teamwork ratings by each team member about others, averaged over time",caption="CDR=Commander, FE=Flight Engineer, MS1/MS3=Mission Specialist, CS=Surgeon, MCC=Mission Control")
ggplot(merged.2,aes(x=Measure,y=sd,group=Rater,color=Rater)) +geom_line(size=1.5) +geom_point(size=3.5) + theme_minimal() +  scale_y_continuous(limits=c(0,5))+
theme(axis.text.x=element_text(size=12)) +
labs(title="Mean variability of teamwork ratings by each team member, averaged over time",
caption="CDR=Commander, FE=Flight Engineer, MS1/MS3=Mission Specialist, CS=Surgeon, MCC=Mission Control")
merged.1 <- merged.byday %>% group_by(measure,rated) %>% summarize(mean=mean(mean,na.rm=T),
sd=mean(sd,na.rm=T))
labs <- c("Trust", "Team\n orientation", "Communications","Leadership","Backup\n behavior","Performance\nmonitoring","Adaptability","Shared\n mental model")
labs <- factor(labs,levels=labs)
merged.1$Measure<- labs[factor(merged.1$measure,levels=c("trust","torient","comm","ls","backup",    "monitor","adapt","smm" ))]
ggplot(merged.1,aes(x=Measure,y=mean,group=rated,color=rated)) +geom_line(size=1.5) +geom_point(size=3.5) + theme_minimal() +  scale_y_continuous(limits=c(0,10))+
theme(axis.text.x=element_text(size=12)) +
labs(title="Teamwork ratings made about each team member by others, averaged over time",caption="CDR=Commander, FE=Flight Engineer, MS1/MS3=Mission Specialist, CS=Surgeon, MCC=Mission Control")
ggplot(merged.1,aes(x=Measure,y=sd,group=rated,color=rated)) +geom_line(size=1.5) +geom_point(size=3.5) + theme_minimal() +  scale_y_continuous(limits=c(0,5))+
theme(axis.text.x=element_text(size=12)) +
labs(title="Mean variability of teamwork ratings by about each team member by others, averaged over time",
caption="CDR=Commander, FE=Flight Engineer, MS1/MS3=Mission Specialist, CS=Surgeon, MCC=Mission Control")
trust.byday <- r.trust.long      %>% group_by(Day,TaskType)  %>% summarize(mean=mean(value,na.rm=T), sd=sd(value,na.rm=T))
torient.byday <- r.torient.long  %>% group_by(Day,TaskType)  %>% summarize(mean=mean(value,na.rm=T), sd=sd(value,na.rm=T))
comm.byday <- r.comm.long        %>% group_by(Day,TaskType)  %>% summarize(mean=mean(value,na.rm=T), sd=sd(value,na.rm=T))
ls.byday <- r.ls.long            %>% group_by(Day,TaskType)  %>% summarize(mean=mean(value,na.rm=T), sd=sd(value,na.rm=T))
backup.byday <- r.backup.long    %>% group_by(Day,TaskType)  %>% summarize(mean=mean(value,na.rm=T), sd=sd(value,na.rm=T))
monitor.byday <- r.monitor.long  %>% group_by(Day,TaskType)  %>% summarize(mean=mean(value,na.rm=T), sd=sd(value,na.rm=T))
adapt.byday <- r.adapt.long      %>% group_by(Day,TaskType)  %>% summarize(mean=mean(value,na.rm=T), sd=sd(value,na.rm=T))
smm.byday <- r.smm.long          %>% group_by(Day,TaskType)  %>% summarize(mean=mean(value,na.rm=T), sd=sd(value,na.rm=T))
all.byday.sd<- data.frame(trust.byday[,1:3],
trust=trust.byday$sd,
torient=torient.byday$sd,
comm=comm.byday$sd,
ls=ls.byday$sd,
backup=backup.byday$sd,
monitor=monitor.byday$sd,
adapt=adapt.byday$sd,
smm=smm.byday$sd
)
all.byday.mean <- data.frame(trust.byday[,1:3],
trust=trust.byday$mean,
torient=torient.byday$mean,
comm=comm.byday$mean,
ls=ls.byday$mean,
backup=backup.byday$mean,
monitor=monitor.byday$mean,
adapt=adapt.byday$mean,
smm=smm.byday$mean
)
par(mfrow=c(1,1))
tmp <- all.byday.mean
matplot(tmp$Day,tmp[,c(4:11)],col=2:9,
type="b",lty=1,pch=16,cex=2.5,
ylab="Mean rating across rating categories",xlab="Mission day",xlim=c(0,120),ylim=c(0,11),
main="Mean ratings over SIRIUS-21 Mission")
matplot(tmp$Day,tmp[,c(4:11)],col="black",
type="b",lty=1,pch=c("T","O","C","L","B","M","A","S"),cex=.85,add=T)
text(tmp$Day,105,tmp$Delay)
grid()
legend(80,5,c("T: Trust","O: Team orientation",
"C: Communication","L: Leadership","B: Backup","M: Performance monitoring",
"A: Adaptability","S: Shared mental model"),
cex=.8,col=2:9,pch=16,pt.cex=1.2)
tmp <- all.byday.sd
matplot(tmp$Day,tmp[,c(4:11)],col=2:9,
type="b",lty=1,pch=16,cex=2.5,
ylab="SD of rating across rating categories",xlab="Mission day",xlim=c(0,120),ylim=c(0,6),
main="SD of ratings over SIRIUS-21 Mission")
matplot(tmp$Day,tmp[,c(4:11)],col="black",
type="b",lty=1,pch=c("T","O","C","L","B","M","A","S"),cex=.85,add=T)
text(tmp$Day,55,tmp$Delay)
grid()
legend(20,6,c("T: Trust","O: Team orientation",
"C: Communication","L: Leadership","B: Backup","M: Performance monitoring",
"A: Adaptability","S: Shared mental model"),
cex=.8,col=2:9,pch=16,pt.cex=1.2)
library(corrplot)
labels <- c("CDR","FE","MS1","MS3","CS","MCC")
r.trust <- select(ratings,Day,TaskType,Rater, starts_with("Trust"))
par(mfrow=c(2,4),xpd=NA)
cc1 <-cor(select(ratings,starts_with("Trust")), use="pairwise.complete",method='spearman')
colnames(cc1) <- rownames(cc1) <- labels
cc2 <-cor(select(ratings,starts_with("Torient")), use="pairwise.complete",method='spearman')
colnames(cc2) <- rownames(cc2) <- labels
cc3 <-cor(select(ratings,starts_with("Comm")), use="pairwise.complete",method='spearman')
colnames(cc3) <- rownames(cc3) <- labels
cc4 <-cor(select(ratings,starts_with("LS")), use="pairwise.complete",method='spearman')
colnames(cc4) <- rownames(cc4) <- labels
cc5 <-cor(select(ratings,starts_with("Backup")), use="pairwise.complete",method='spearman')
colnames(cc5) <- rownames(cc5) <- labels
cc6 <-cor(select(ratings,starts_with("Monitor")), use="pairwise.complete",method='spearman')
colnames(cc6) <- rownames(cc6) <- labels
cc7 <-cor(select(ratings,starts_with("Adapt")), use="pairwise.complete",method='spearman')
colnames(cc7) <- rownames(cc7) <- labels
cc8 <-cor(select(ratings,starts_with("SMM")), use="pairwise.complete",method='spearman')
colnames(cc8) <- rownames(cc8) <- labels
corrplot(cc1,title =paste("SIRIUS 21 Mission","Trust"), order="orig",mar=c(0,0,5,1),type="upper",is.corr=T,diag=F)
corrplot(cc2,title="Team orientation",order="orig",mar=c(0,0,5,1),type="upper",is.corr=T,diag=F)
corrplot(cc3,main="Communication",order="orig",mar=c(0,0,5,1),type="upper",is.corr=T,diag=F)
corrplot(cc4,main="Leadership",order="orig",mar=c(0,0,5,1),type="upper",is.corr=T,diag=F)
corrplot(cc5,main="Backup",order="orig",mar=c(0,0,5,1),type="upper",is.corr=T,diag=F)
corrplot(cc6,main="Performance monitoring",order="orig",mar=c(0,0,5,1),type="upper",is.corr=T,diag=F)
corrplot(cc7,main="Adaptability",order="orig",mar=c(0,0,5,1),type="upper",is.corr=T,diag=F)
corrplot(cc8,main="Shared MM",order="orig",mar=c(0,0,5,1),type="upper",is.corr=T,diag=F)
cat("\n\n\\pagebreak\n")
par(mfrow=c(1,1))
avg <-(cc1+cc2+cc3+cc4+cc5+cc6+cc7+cc8)/8
rownames(avg) <- colnames(avg) <- labels
par(mar=c(2,5,20,2))
corrplot(avg, main=paste("Overall similarity of ratings about each role\n(SIRIUS 21 Mission)",sep=""),order="orig",mar=c(0,0,5,1),
is.corr=T,diag=F,type="upper")
cat("\n\n\\pagebreak\n")
library(corrplot)
library(reshape2)
labels <- c("CDR","FE","MS1","MS2","MS3","CS","MCC")
##these are the ratings provided by each rater/role about each rater.  Both include MCC.
r1 <-select(ratings,Day,Rater,starts_with("Trust"))
r2 <-select(ratings,Day,Rater,starts_with("Torient"))
r3 <-select(ratings,Day,Rater,starts_with("Comm"))
r4 <-select(ratings,Day,Rater,starts_with("LS"))
r5 <-select(ratings,Day,Rater,starts_with("Backup"))
r6 <-select(ratings,Day,Rater,starts_with("Monitor"))
r7 <-select(ratings,Day,Rater,starts_with("Adapt"))
r8 <-select(ratings,Day,Rater,starts_with("SMM"))
colnames(r1)<- c("Day","Role",labels)
colnames(r2)<- c("Day","Role",labels)
colnames(r3)<- c("Day","Role",labels)
colnames(r4)<- c("Day","Role",labels)
colnames(r5)<- c("Day","Role",labels)
colnames(r6)<- c("Day","Role",labels)
colnames(r7)<- c("Day","Role",labels)
colnames(r8)<- c("Day","Role",labels)
##now, transform these to
r1.t <- melt(r1,id.var=c("Day","Role")) %>% group_by(Role) %>% dcast(Day+variable~Role)
r2.t <- melt(r2,id.var=c("Day","Role")) %>% group_by(Role) %>% dcast(Day+variable~Role)
r3.t <- melt(r3,id.var=c("Day","Role")) %>% group_by(Role) %>% dcast(Day+variable~Role)
r4.t <- melt(r4,id.var=c("Day","Role")) %>% group_by(Role) %>% dcast(Day+variable~Role)
r5.t <- melt(r5,id.var=c("Day","Role")) %>% group_by(Role) %>% dcast(Day+variable~Role)
r6.t <- melt(r6,id.var=c("Day","Role")) %>% group_by(Role) %>% dcast(Day+variable~Role)
r7.t <- melt(r7,id.var=c("Day","Role")) %>% group_by(Role) %>% dcast(Day+variable~Role)
r8.t <- melt(r8,id.var=c("Day","Role")) %>% group_by(Role) %>% dcast(Day+variable~Role)
cc1.t <- cor(r1.t[,3:7],use="pairwise.complete",method="spearman")
cc2.t <- cor(r2.t[,3:7],use="pairwise.complete",method="spearman")
cc3.t <- cor(r3.t[,3:7],use="pairwise.complete",method="spearman")
cc4.t <- cor(r4.t[,3:7],use="pairwise.complete",method="spearman")
cc5.t <- cor(r5.t[,3:7],use="pairwise.complete",method="spearman")
cc6.t <- cor(r6.t[,3:7],use="pairwise.complete",method="spearman")
cc7.t <- cor(r7.t[,3:7],use="pairwise.complete",method="spearman")
cc8.t <- cor(r8.t[,3:7],use="pairwise.complete",method="spearman")
if(1)
{
par(mfrow=c(2,4),xpd=F)
corrplot(cc1.t,title =paste("SIRIUS-21 Mission","Trust"),order="orig",mar=c(0,0,1,1),type="upper",is.corr=T,diag=F)
corrplot(cc2.t,title="Team orientation",order="orig",mar=c(0,0,1,1),type="upper",is.corr=T,diag=F)
corrplot(cc3.t,main="Communication",order="orig",mar=c(0,0,1,1),type="upper",is.corr=T,diag=F)
corrplot(cc4.t,main="Leadership",order="orig",mar=c(0,0,1,1),type="upper",is.corr=T,diag=F)
corrplot(cc5.t,main="Backup",order="orig",mar=c(0,0,1,1),type="upper",is.corr=T,diag=F)
corrplot(cc6.t,main="Performance monitoring",order="orig",mar=c(0,0,1,1),type="upper",is.corr=T,diag=F)
corrplot(cc7.t,main="Adaptability",order="orig",mar=c(0,0,1,1),type="upper",is.corr=T,diag=F)
corrplot(cc8.t,main="Shared MM",order="orig",mar=c(0,0,1,1),type="upper",is.corr=T,diag=F)
cat("\n\n\\pagebreak\n")
}
avg <-( cc1.t+cc2.t+cc3.t+cc4.t+cc5.t+cc6.t+cc7.t+cc8.t)/8
par(mfrow=c(1,1))
corrplot(avg, main=paste("Correlation of ratings from each rater\n (SIRIUS Mission)",sep=""),order="orig",mar=c(1,0,2,1),
is.corr=T,diag=T,type="upper")
library(pwr)
install.packages("pwr")
library(pwr)
#Test Type A:
pab_saA <- t.test(TeamA_sa$TotalA,TeamB_sa$TotalA,alternative="greater",paired=F)
library(dplyr)
?filter
x <- 1:100
plot(x,(100-x)/2)
library(elo)
#library(EloRating) #Animal dominance
#library(EloChoice) # Visual stimuli parwise comparisons
#library(EloOptimized)
data(tournament)
str(tournament)
elo.run(score(points.Home, points.Visitor) ~ team.Home + team.Visitor, data = tournament, k = 20)
x <- elo.run(score(points.Home, points.Visitor) ~ team.Home + team.Visitor +
k(20*log(abs(points.Home - points.Visitor) + 1)), data = tournament)
logodds <- function(p) {log(p/(1-p))}
logit <- function(x) {(exp(x)/(1 + exp(x)))}
##this should work for an entire vectorized list, so p1 and p2 can be a series of
##opponents over which skill is static
play <- function(p1,p2)
{
skilldelta <- t(t(skill[p1,,drop=F] - skill[p2,,drop=F]) * skillWeights)
##make winner depend on a weighted average
sumdelta <- rowSums(skilldelta)
advantage <- sumdelta #could use other strategies for multidimensional things.
winprob <- logit(advantage)
wins <- runif(length(winprob))<winprob
return (wins)
}
MatchMake<- function(p1)
{
pOthers <- sample(1:nPlayers,nAvailable)
pOthers <- pOthers[pOthers!=p1]
elo1 <-Elo[p1]
elo2 <- Elo[pOthers]
rankSim <- pOthers[order((elo1 - elo2)^2)]
prob <- 1/(1:length(pOthers))
prob <- prob/sum(prob)
pairing <- sample(rankSim,1,prob=prob)
pairing
}
updateElo <- function(elo1, elo2, p1Win,k=32)
{
change <- (elo2-elo1)
expectedWin <- 1/(1+10^(change/400))
diff <- k * (p1Win-expectedWin)
c(elo1+diff,elo2-diff)
}
nPlayers <- 250
nSkills <- 10
nAvailable <- 50   ##how many are available at any given time. This will actually
## be a distribution that is likely to be related to how good you are,
## how much fun you are having, etc.
##normal:
skill <- matrix(10+rnorm(nPlayers*nSkills),nrow=nPlayers)
##skewed
#skill <- matrix(exp(2+rnorm(nPlayers*nSkills)),nrow=nPlayers)
##Let's rank-order them by weighted skill.
skillWeights <- rep(0,nSkills)
skillWeights[1] <- 1
skill <- skill[order(-skill %*% skillWeights),]
##starting ELO of 1000
Elo <- 999.5+runif(nPlayers)
sampleDist <-rep(1,nPlayers)/nPlayers
sampleDist <- (nPlayers:1)
games <- 40000
eloHist <- matrix(0,nrow=games,ncol=nPlayers)
for(i in 1:games)
{
p1 <- sample(1:nPlayers,1,prob=sampleDist)
p2 <- MatchMake(p1)
eloPre <- Elo[c(p1,p2)]
outcome <- play(p1,p2)
eloPost <- updateElo(Elo[p1],Elo[p2],outcome)
Elo[c(p1,p2)] <- eloPost
#print(c(p1,p2,outcome,eloPre,eloPost))
eloHist[i,] <- Elo
}
matplot(eloHist,type="l",col="black",lty=1)
plot(skill[,1],Elo)
hist(Elo,breaks=25,xlim=c(0,2000))
skill <- matrix(exp(2+runif(nPlayers*nSkills)),nrow=nPlayers)
nPlayers <- 250
nSkills <- 10
nAvailable <- 50   ##how many are available at any given time. This will actually
## be a distribution that is likely to be related to how good you are,
## how much fun you are having, etc.
##normal:
skill <- matrix(10+rnorm(nPlayers*nSkills),nrow=nPlayers)
##skewed
#skill <- matrix(exp(2+rnorm(nPlayers*nSkills)),nrow=nPlayers)
skill <- matrix(exp(2+runif(nPlayers*nSkills)),nrow=nPlayers)
##Let's rank-order them by weighted skill.
skillWeights <- rep(0,nSkills)
skillWeights[1] <- 1
skill <- skill[order(-skill %*% skillWeights),]
##starting ELO of 1000
Elo <- 999.5+runif(nPlayers)
sampleDist <-rep(1,nPlayers)/nPlayers
sampleDist <- (nPlayers:1)
games <- 40000
eloHist <- matrix(0,nrow=games,ncol=nPlayers)
for(i in 1:games)
{
p1 <- sample(1:nPlayers,1,prob=sampleDist)
p2 <- MatchMake(p1)
eloPre <- Elo[c(p1,p2)]
outcome <- play(p1,p2)
eloPost <- updateElo(Elo[p1],Elo[p2],outcome)
Elo[c(p1,p2)] <- eloPost
#print(c(p1,p2,outcome,eloPre,eloPost))
eloHist[i,] <- Elo
}
matplot(eloHist,type="l",col="black",lty=1)
plot(skill[,1],Elo)
hist(Elo,breaks=25,xlim=c(0,2000))
install.packages("rmarkdown")
plot(x,x)
?emmeans
library(emmeans)
??eammeans
?emmeans
20000000/200000
4000000/17
set.seed(100)
c(76834,-8162)/8
knitr::opts_chunk$set(echo = TRUE)
cardat <- read.table(text="age gender  type origin origin.last carval carval.last
1   34      M Truck     US      Europe  16400       15900
2   31      M Truck Europe      Europe  16900       15800
3   47      F Sedan Europe      Europe  18800       16600
4   21      F   SUV  Japan       Japan  16000       16200
5   42      F Truck     US          US  16800       16500
6   43      F Truck  Japan       Japan  17200       16200
7   60      M Truck  Japan       Japan  19900       20800
8   37      F Sedan  Japan       Japan  17100       16200
9   46      M   SUV     US          US  16900       16700
10  27      F Sedan     US       Japan  16200       15600
11  50      M   SUV     US      Europe  18800       17800
12  64      F Sedan  Japan       Japan  50700       18000
13  33      M   SUV     US          US  16500       15700
14  39      M   SUV     US          US  17000       17600
15  58      F Truck     US          US  19400       16600
16  53      M Sedan  Japan          US  19200       17700
17  29      F Sedan  Japan          US  16300       15600
18  37      M Truck  Japan      Europe  17300       16200
19  37      M Truck     US          US  18200       16600
20  54      M Sedan     US      Europe  24500       16400
21  46      M Truck  Japan      Europe  18000       17700
22  55      F Truck     US          US  28900       17900
23  46      F Sedan  Japan       Japan  16600       16200
24  57      F Truck Europe      Europe  24300       16600
25  40      F Truck     US          US  16800       16700
26  27      F Sedan  Japan          US  16900       15600
27  58      M Sedan     US          US  20300       19000
28  64      M Sedan     US          US  40600       23600
29  47      F Sedan     US          US  18400       19800
30  32      F   SUV     US      Europe  15900       15800
31  43      F Sedan  Japan      Europe  17200       16300
32  66      M   SUV     US          US  19100       30200
33  36      F   SUV  Japan      Europe  16900       15800
34  68      F Sedan  Japan      Europe  69300       18500
35  54      M Sedan     US          US  17000       17200
36  64      F   SUV     US          US  34900       42900
37  27      F Truck     US       Japan  15800       15500
38  51      F Sedan  Japan       Japan  29000       17900
39  69      F Sedan Europe       Japan  54400       18700
40  25      M Truck  Japan       Japan  15800       15400",header=T)
chisq.test(table(cardat$origin,cardat$origin.last))
setwd("~/Dropbox/courses/5210-2022/problemsets/ps6")
ns
ns()
library(gam)
ns
dat <- read.csv("ptracker.csv")[200:500,]
setwd("~/Dropbox/courses/5210-2022/web-5210/psy5210/Projects/Chapter13")
dat <- read.csv("ptracker.csv")[200:500,]
time <- 117:168*100  ##Regularize in tenths of seconds
mx <- loess(posx~time,span=.1,data=dat)
my <- loess(posy~time,span=.1,data=dat)
par(mfrow=c(1,3))
plot(dat$time,dat$posx,main="X position")
points(time,predict(mx,data.frame(time=time)),type="l")
plot(dat$time,dat$posy,main="Y position")
points(time,predict(my,data.frame(time=time)),type="l")
plot(dat$posx,dat$posy,type="b",main="X and Y")
points(predict(mx,data.frame(time=time)),
predict(my,data.frame(time=time) ),
col="red",type="o",pch=1,cex=.3)
cw <- ChickWeight
cw
cw[1:5,]
g <-  gam(Diet~s(time)+ Diet)
cw <- ChickWeight
g <-  gam(Diet~s(time)+ Diet)
g <-  gam(weight~s(time:Chick)+ Diet,data=cw)
g <-  gam(weight~s(time)*Chick+ Diet,data=cw)
library(gam)
cw <- ChickWeight
g <-  gam(weight~s(Time)*Chick+ Diet,data=cw)
hg
g
summary(g)
plot(g)
g$fitted.values\
g$fitted.values
ggplot(cw,aes(x=Day,y=fit,group=Diet,color=Diet)) +geom_point() + geom_line()
library(ggplot2)
ggplot(cw,aes(x=Day,y=fit,group=Diet,color=Diet)) +geom_point() +
theme_bw())
library(ggplot2)
ggplot(cw,aes(x=Day,y=fit,group=Diet,color=Diet)) +geom_point() + theme_bw()
cw[1:10,]
cw$fit <- g$fitted.values
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit,group=Diet,color=Diet)) +geom_point() + theme_bw()
summary(g)
library(gam)
cw <- ChickWeight
g <-  gam(weight~s(Time)+ Diet,data=cw)
cw$fit <- g$fitted.values
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit,group=Diet,color=Diet)) +geom_point() + theme_bw()
plot(g)
mx <- lm(posx~m(time),data=dat)
mx <- lm(posx~s(time),data=dat)
mx
plot(mx)
mx <- lm(posx~ns(time,df=5),data=dat)
plot(mx)
plot(dat$time,mx$fit)
plot(dat$time,dat$posx,type="l")
points(dat$time,mx$fit,col="red")
mx5 <- lm(posx~ns(time,df=5),data=dat)
mx10 <- lm(posx~ns(time,df=10),data=dat)
mx15 <- lm(posx~ns(time,df=15),data=dat)
plot(dat$time,dat$posx,type="l")
points(dat$time,mx5$fit,col="red")
points(dat$time,mx10$fit,col="gold")
points(dat$time,mx15$fit,col="navy")
par(mfrow=c(1,1))
plot(dat$time,dat$posx,type="l")
points(dat$time,mx5$fit,col="red")
points(dat$time,mx10$fit,col="gold")
points(dat$time,mx15$fit,col="navy")
dat[1:5,]
?ns
g <-  gam(weight~s(Time)*Chick+ Diet,data=cw)
cw$fit2 <- g$fitted.values
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit2,group=Diet,color=Diet)) +geom_point() + theme_bw()
summary(g)
plot(g)
ggplot(cw,aes(x=Time,y=fit2,group=Diet,color=Diet)) +geom_point() + theme_bw()
ggplot(cw,aes(x=Time,y=fit2,group=Diet,color=Diet)) +geom_point() + geom_line()+ theme_bw()
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit2,group=Chick,color=Diet)) +geom_point() + geom_line()+ theme_bw()
summary(g)
#plot(g)
summary(g)
g <-  gam(weight~s(Time)+ Diet,data=cw)
cw$fit <- g$fitted.values
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit,group=Diet,color=Diet,group=Chick)) +geom_point() +geom_line() theme_bw()
library(gam)
cw <- ChickWeight
g <-  gam(weight~s(Time)+ Diet,data=cw)
cw$fit <- g$fitted.values
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit,group=Diet,color=Diet,group=Chick)) +geom_point() +geom_line() +theme_bw()
summary(g)
plot(g)
cw <- ChickWeight
g <-  gam(weight~s(Time)+ Diet,data=cw)
cw$fit <- g$fitted.values
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit,group=Diet,color=Diet,group=Chick)) +geom_point() +geom_line() +theme_bw()
summary(g)
plot(g)
g
plot(g)
?gam
bs(cw$weight,cw$Time,df=10)
plot(bs(cw$weight,cw$Time,df=10))
plot(bs(cw$weight,cw$Time,df=10)[1])
require(stats); require(graphics)
bs(women$height, df = 5)
summary(fm1 <- lm(weight ~ bs(height, df = 5), data = women))
## example of safe prediction
plot(women, xlab = "Height (in)", ylab = "Weight (lb)")
ht <- seq(57, 73, length.out = 200)
lines(ht, predict(fm1, data.frame(height = ht)))
?ns
data(gam.data)
Gam.object <- gam(y~x+z, data=gam.data)
step.object <-step.Gam(Gam.object, scope=list("x"=~1+x+s(x,4)+s(x,6)+s(x,12),"z"=~1+z+s(z,4)))
## Not run:
step.Gam(Gam.object, scope=list("x"=~1+x+s(x,4)+s(x,6)+s(x,12),"z"=~1+z+s(z,4)),parallel=TRUE)
best <- step.Gam(Gam.object, scope=list("x"=~1+x+s(x,4)+s(x,6)+s(x,12),"z"=~1+z+s(z,4)))
best
plot(best)
?gam
summary(best)
library(gam)
cw <- ChickWeight
g <-  gam(weight~s(Time)+ Diet,data=cw)
cw$fit <- g$fitted.values
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit,group=Diet,color=Diet,group=Chick)) +geom_point() +geom_line() +theme_bw()
summary(g)
plot(g)
plot(g)
ggplot(cw,aes(x=Time,y=fit,group=Diet,color=Diet,group=Chick)) +geom_point() +geom_line() +theme_bw()
g <-  gam(weight~lo(Time)*Chick+ Diet,data=cw)
cw$fit2 <- g$fitted.values
library(ggplot2)
ggplot(cw,aes(x=Time,y=fit2,group=Chick,color=Diet)) +geom_point() + geom_line()+ theme_bw()
summary(g)
g <-  lm(weight~lo(Time)*Chick+ Diet,data=cw)
summary(g)
summary(g)
g
g <-  lm(weight~ns(Time,10)*Chick+ Diet,data=cw)
summary(g)
g <-  lm(weight~ns(Time,5)++ Diet,data=cw)
g <-  lm(weight~ns(Time,5)+ Diet,data=cw)
g
anova(g)
summary(g)