tyuma.df<-read.table("exam1.tyuma.data",header=T)
##
cat("\n\n")
cat("Trace output for unweighted fit of Tyuma data to MWC model\n")
ty.mwc<-nls(y ~
(L/Kt*p*(1+p/Kt)^3+p/Kr*(1+p/Kr)^3)/(L*(1+p/Kt)^4+(1+p/Kr)^4),tyuma.df,list(L=3.e5,Kt=10,Kr=.1),control=nls.control(maxiter=100,minFactor=.001),trace=T)
##
cat("\n\n")
cat("Trace output for weighted fit of Tyuma data to MWC model\n")
ty.mwc.wt<-nls(rep(0,length(y)) ~ wt.mwc(y, p, L, Kt,
Kr),tyuma.df,list(L=3.e5,Kt=10,Kr=.1),control=nls.control(maxiter=100,minFactor=.001),trace=T)
##
## unweighted fit to adair equation dumps -- must weight the residuals
cat("\n\n")
cat("Trace output for weighted fit of Tyuma data to Adair model\n")
ty.ad.wt<-nls(rep(0,length(y)) ~ wt.adair(y, p, K1,K2,K3,K4),tyuma.df,list(K1=.024,K2=.0743,K3=.0858,K4=7.37),control=nls.control(maxiter=100,minFactor=.001),trace=T)
##
## attach dataframe so can address input data as "p" and "y"
attach(tyuma.df)
##
## construct variables holding calculated values of y
y.mwc<-fitted.values(ty.mwc)
y.mwc.wt<-(y - residuals(ty.mwc.wt)*y*(1-y))
y.ad.wt<-(y - residuals(ty.ad.wt)*y*(1-y))
sumdelsq.mwc<-sum((y - y.mwc)^2)
sumdelsq.mwc.wt<-sum((y - y.mwc.wt)^2)
sumdelsq.ad.wt<-sum((y - y.ad.wt)^2)
sumdelsq.vec<-c(mwc=sumdelsq.mwc,mwc.wt=sumdelsq.mwc.wt,ad.wt=sumdelsq.ad.wt)
sumdelsq.wt.vec<-c(mwc=sum(((y - y.mwc)/(y*(1-y)))^2), mwc.wt=sum(((y -
y.mwc.wt)/(y*(1-y)))^2), ad.wt=sum(((y - y.ad.wt)/(y*(1-y)))^2))
cat("\n\n")
cat("Summary for unweighted fit of Tyuma data to MWC model\n")
print(summary(ty.mwc))
cat("\n\n")
cat("Summary for weighted fit of Tyuma data to MWC model\n")
print(summary(ty.mwc.wt))
cat("\n\n")
cat("Summary for weighted fit of Tyuma data to Adair model\n")
print(summary(ty.ad.wt))
cat("\n\n")
cat("\n\nUnweighted variance = sum((y - ycalc)^2)/df, \nfor unweighted and weighted fits\n\n")
print(sumdelsq.vec/c(17,17,16))
cat("\n")
cat("\n\nWeighted variance = sum(((y - ycalc)/(y*(1-y)))^2)/df, \nfor unweighted and weighted fits\n\n")
print(sumdelsq.wt.vec/c(17,17,16))
##
plot(p,y,pch=3)
points(p,y.mwc,pch=22,col=2)
points(p,y.mwc.wt,pch=24,col=3)
lines(p,y.mwc,lty=3,col=2)
lines(p,y.mwc.wt,lty=1,col=3)
#lines(smooth.spline(p,y.mwc,spar=.1),lty=3)
#lines(smooth.spline(p,y.mwc.wt,spar=.1),lty=1)
#lines(lowess(p,y.mwc,f=.1),lty=3)
#lines(lowess(p,y.mwc.wt,f=.1),lty=1)
title("Data of Tyuma et al.\nFractional satn, y, vs. Oxygen pressure, p")
text(25,.4,"Fits to MWC equation:\n\ndashed line, uniformly weighted\n\nsolid line, weighted by 1/(y*(1 - y))")
cat("type <return> to continue: ")
readline()
##
plot(p,y,pch=3,log="xy")
points(p,y.mwc,pch=22,col=2)
points(p,y.mwc.wt,pch=24,col=3)
lines(p,y.mwc,lty=3,col=2)
lines(p,y.mwc.wt,lty=1,col=3)
#lines(smooth.spline(p,y.mwc,spar=.4),lty=3)
#lines(smooth.spline(p,y.mwc.wt,spar=.4),lty=1)
#lines(lowess(p,y.mwc),lty=3)
#lines(lowess(p,y.mwc.wt),lty=1)
title("Data of Tyuma et al.\nlog(y)) vs. log(p)")
text(10,.05,"Fits to MWC equation:\n\ndashed line, uniformly weighted\n\nsolid line, weighted by 1/(y*(1 - y))")
cat("type <return> to continue: ")
readline()
##
plot(p,y/(1-y),pch=3,log="xy")
points(p,y.mwc/(1-y.mwc),pch=22,col=2)
points(p,y.mwc.wt/(1-y.mwc.wt),pch=24,col=3)
lines(p,y.mwc/(1-y.mwc),lty=3,col=2)
lines(p,y.mwc.wt/(1-y.mwc.wt),lty=1,col=3)
#lines(smooth.spline(p,y.mwc/(1-y.mwc),spar=.4),lty=3)
#lines(smooth.spline(p,y.mwc.wt/(1-y.mwc.wt),spar=.4),lty=1)
#lines(lowess(p,y.mwc/(1-y.mwc),f=.5),lty=3)
#lines(lowess(p,y.mwc.wt/(1-y.mwc.wt),f=.5),lty=1)
title("Data of Tyuma et al.\nlog(y/(1 - y)) vs. log(p)")
text(10,.1,"Fits to MWC equation:\n\ndashed line, uniformly weighted\n\nsolid line, weighted by 1/(y*(1 - y))")
cat("type <return> to continue: ")
readline()
##
plot(p,y/(1-y),pch=3,log="xy")
points(p,y.ad.wt/(1-y.ad.wt),pch=22,col=2)
points(p,y.mwc.wt/(1-y.mwc.wt),pch=24,col=3)
lines(p,y.ad.wt/(1-y.ad.wt),lty=3,col=2)
lines(p,y.mwc.wt/(1-y.mwc.wt),lty=1,col=3)
#lines(smooth.spline(p,y.ad.wt/(1-y.ad.wt),spar=.4),lty=3)
#lines(smooth.spline(p,y.mwc.wt/(1-y.mwc.wt),spar=.4),lty=1)
#lines(lowess(p,y.ad.wt/(1-y.ad.wt),f=.5),lty=3)
#lines(lowess(p,y.mwc.wt/(1-y.mwc.wt),f=.5),lty=1)
title("Data of Tyuma et al.\nlog(y/(1 - y)) vs. log(p)")
text(10,.1,"Fits weighted by 1/(y*(1 - y)):\n\ndashed line, Adair eqn\n\nsolid line, MWC eqn")
##
## xgobi(cbind(tyuma.df,y.mwc=y.mwc, y.mwc.wt=y.mwc.wt, y.adair.wt=y.ad.wt))
##