! -------------------------------------
! Author: David W. Gerbing
!         Portland State University
!         davidg@sba.pdx.edu
! Date:   Dec 30, 2010
! Purpose: cgi program for R regression
! -------------------------------------
program regression
            
implicit none

integer :: i, j, nchars, nchars2, p, nX, nsubs=0
logical :: infl_flag, pred_flag, cor_flag, best_flag, resid_flag
character (len=1) :: c, c1, c2, c3, c4, c5
character (len=15) :: X(9), Y, wherestr, dummy
character (len=30) :: Ystr, Xstr(9), cp
character (len=100) :: datastr, txt1, txt2, ctitle
character (len=200) :: line, line2, fname
character (len=400) :: cm


! -------------------------------------------------------------
! ------------------ Read and Parse CGI Data ------------------
! -------------------------------------------------------------
read(*,'(a)') line
! write(*,'(a)') line

! see if data file is local or on the net
p = index(line,"=")
cp = line(p+1:p+1)
if (cp == "l") then
  wherestr = "local"
else
	wherestr = "net"
end if

! read number of X vars and convert to number
p = index(line,"nX")
cp = line(p+3:p+3)
read(cp,'(i3)') nX

! prepare line for fortran csv read that parses line
! also replace the browser add + for a space with a .
nchars = len(trim(line))
do i = 1, nchars
  if (line(i:i) == '&') line(i:i) = ','
  if (line(i:i) == '+') line(i:i) = '.'
end do

!write(*,'(a)') line
!do i = 1, nchars
!  write(*,'(i3,1x,a)') i, line(i:i)
!end do

! browser escape character processing
do i = 1, nchars
	c = line(i:i)
  if (c == "%") then
    if(line(i+1:i+2) == "3A") line(i:i) = ':'
    if(line(i+1:i+2) == "2F") line(i:i) = '@'  ! add / later
    if(line(i+1:i+2) == "7E") line(i:i) = '~'
  end if
  if (c .ne. line(i:i)) then  ! substitution made
		nsubs = nsubs + 1
		line(i+1:i+2) = "  "
	end if
end do

! remove blanks and then remove any extra chars at end of line
j = 0
do i=1,nchars
  c=line(i:i)
	if (c /= ' ') then
		j=j+1
		line2(j:j) = c
  end if
end do
nchars2 = nchars - nsubs*2
do i=nchars2+1, len(line2)
  line2(i:i) = ' '
end do

! ----------------------------------------------------------
! line2 is the parsed, processed character string

! check for checked checkboxes
if (index(line2,"resid") == 0) then
	resid_flag = .false.
else
	resid_flag = .true.
end if
if (index(line2,"pred") == 0) then
	pred_flag = .false.
else
	pred_flag = .true.
end if
if (index(line2,"cor") == 0) then
	cor_flag = .false.
else
	cor_flag = .true.
end if
if (index(line2,"best") == 0) then
	best_flag = .false.
else
	best_flag = .true.
end if

! write(*,'(a)') line2
!do i = 1, nchars2
! write(*,'(i3,1x,a)') i, line2(i:i)
!end do

! read data file name if applicable, then Y and X vars
if (wherestr == "local") then
	read(unit=line2, fmt=*) dummy, Ystr, (Xstr(i), i=1,nX)
	datastr = ""
else
	read(unit=line2, fmt=*) dummy, datastr, Ystr, (Xstr(i), i=1,nX)
	do i = 1, len(trim(datastr))
		if (datastr(i:i) == '@') datastr(i:i) = '/'
	end do
end if
! write(*,'(a)') line2

! write(*,*) '*** ', datastr, Ystr, (Xstr(i), i=1,nX)
! parse each name-value pair into just the value
if (wherestr == "local") then
	fname = ""
else
	fname = datastr(index(datastr,"=")+1:len(trim(datastr)))
end if
Y = Ystr(index(Ystr,"=")+1:len(trim(Ystr)))
do i=1,nX
	X(i) = Xstr(i)(index(Xstr(i),"=")+1:len(trim(Xstr(i))))
enddo


! -------------------------------------------------------------
! ------------- Construct R Regression Program ----------------
! -------------------------------------------------------------

! num to char, used to designate columns
write(c1,'(i1)') nX+2; write(c2,'(i1)') nX+3; write(c3,'(i1)') nX+4
write(c4,'(i1)') nX+5; write(c5,'(i1)') nX+6


! must precede all other write statements to stdout
write(*,'(a,//)') "Content-type:text/html"
! write(*,'(a,/)') line
! write(*,'(a)') line2
write(*,'(a)') '<html>'
write(*,'(a)') '<body>'

write(*,'(a)') '<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"'
write(*,'(a)') '       "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">'
write(*,'(a)') '<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en"><head>'
write(*,'(a)') '<title>R Regression</title>'
write(*,'(a)') '<meta http-equiv="Content-type" content="text/html; charset=iso-8859-1" />'
write(*,'(a)') '<link rel="stylesheet" href="regR.css" type="text/css" />'
write(*,'(a)') '</head>'
write(*,'(a)') '<body>'

write(*,'(a)') 
write(*,'(a)') '<div class="ltgr">'
write(*,'(a)') '# *-----------------------------------*<br />'
write(*,'(a)') '# R program for OLS regression analysis<br />'
write(*,'(a)') '# *-----------------------------------*<br />'
write(*,'(a)') '# option - encourage R to avoid scientific notation<br />'
write(*,'(a)') '</div>'
write(*,'(a)') '<span class="b">options</span>(scipen=10)<br />'
write(*,'(a)') '<br />'

if (nX == 1 .AND. best_flag) then
	write(*,'(a)') '<div class="ltgr">'
	write(*,'(a)') '# Best subset analysis not applicable when only one predictor.<br /><br />'
	write(*,'(a)') '</div>'
	best_flag = .FALSE.
end if

write(*,'(a)') '<div class="ltgr">'
write(*,'(a)') '# ------------------------<br />'
write(*,'(a)') '# Read and attach the data<br />'
write(*,'(a)') '# ------------------------<br />'
write(*,'(a)') '# data needs to be a '
write(*,'(a)') '<a href="http://web.pdx.edu/~gerbing/Rweb/csv.pdf">.csv<a> text file.<br />'
write(*,'(a)') '</div>'

if(fname .ne. '') then
	write(*,'(a)') 'rfile = "<span class="r">'//trim(fname)//'</span>"<br />'
	write(*,'(a)') 'mydata &lt;- <span class="b">read.csv</span>(file=rfile)<br />'
else
	write(*,'(a)') 'mydata &lt;- <span class="b">read.csv</span>(file.choose())<br />'
endif
write(*,'(a)') '<span class="b">attach</span>(mydata)<br />'
write(*,'(a)') '<br />'
	write(*,'(a)') '<div class="ltgr">'
write(*,'(a)') '# -----------------------------------------------------<br />'
write(*,'(a)') '# List Variable Names and first and last 3 rows of data<br />'
write(*,'(a)') '# ----------------------------------------------------<br />'
write(*,'(a)') '# variable names are what R read from first row of the data file<br />'
write(*,'(a)') '# R replaces each blank in a name with a period<br />'
write(*,'(a)') '</div>'
write(*,'(a)') '<span class="b">head</span>(mydata, n=3)<br />'
write(*,'(a)') '<span class="b">tail</span>(mydata, n=3)<br />'
write(*,'(a)') '<br />'

if (nX == 1) then  ! Need to sort for scatterplot before doing reg analysis
	write(*,'(a)') '<div class="ltgr">'
	write(*,'(a)') '#  Sort according to predictor variable for prediction interval plot<br />'
	write(*,'(a)') '</div>'
	write(*,'(a)') '<span class="b">detach</span>(mydata)<br />'
	txt1 = 'o <- <span class="b">order</span>(mydata[,"<span class="r">'
	write(*,'(a)') trim(txt1)//trim(X(1))//'</span>"], decreasing=FALSE)<br />'
	write(*,'(a)') 'mydata <- mydata[o,]<br />'
	write(*,'(a)') '<span class="b">attach</span>(mydata)<br />'
	write(*,'(a)') '<br />'
end if

write(*,'(a)') 
write(*,'(a)') '<div class="ltgr">'
write(*,'(a)') '# ------------------------------------------------------<br />'
write(*,'(a)') '# Linear Model function - do the OLS regression analysis<br />'
write(*,'(a)') '# ------------------------------------------------------<br />'
write(*,'(a)') '# lm function is of the general form lm(Y ~ X1 + X2 + X3)<br />'
write(*,'(a)') '# Y is response variable, X is one or more predictors<br />'
write(*,'(a)') '# store results in the object called lm.out<br />'
write(*,'(a)') '</div>'
cm = 'lm.out &lt;- <span class="b">lm</span>(<span class="r">'//trim(Y)//'</span>'
cm = trim(cm) // ' ~ <span class="r">'//trim(X(1))//'</span>'
if (nX > 1) then
	do i=2,nX
		cm = trim(cm) //  ' + <span class="r">'//trim(X(i))//'</span>'
	end do
end if
write(*,'(a)') trim(cm)// ')<br /><br />'

write(*,'(a)') 
write(*,'(a)') '<div class="ltgr">'
write(*,'(a)') '# ----------------------------------------<br />'
write(*,'(a)') '# Basic Results of the regression analysis<br />'
write(*,'(a)') '# ----------------------------------------<br />'
write(*,'(a)') '# summary - estimated model<br />'
write(*,'(a)') '# confint - confidence interval for each estimated coefficient<br />'
write(*,'(a)') '# anova - analysis of variance<br />'
write(*,'(a)') '</div>'

write(*,'(a)') '<span class="b">summary</span>(lm.out)<br />'
write(*,'(a)') '<span class="b">confint</span>(lm.out)<br />'
write(*,'(a)') '<span class="b">anova</span>(lm.out)<br />'
write(*,'(a)') '<br />'

if (cor_flag) then  ! correlations and scatteplot matrix
	write(*,'(a)') 
	write(*,'(a)') '<div class="ltgr">'
	write(*,'(a)') '# ---------------------------------<br />'
	write(*,'(a)') '# Correlations of Y and X variables<br />'
	write(*,'(a)') '# ---------------------------------<br />'
	write(*,'(a)') '</div>'
	cm = '(<span class="r">' // trim(Y) // '</span>,<span class="r">'// trim(X(1)) // '</span>'
	if (nX > 1) then
		do i=2,nX
			cm = trim(cm) // ',<span class="r">' // trim(X(i)) // '</span>'
		end do
	end if
	cm = trim(cm) // ')'
	write(*,'(a)') '<span class="b">print</span>( '
	write(*,'(a)') '<span class="b">cor</span>(<span class="b">cbind</span>'//trim(cm)//'), '
	write(*,'(a)') 'digits=2)<br /><br />'
end if

if (resid_flag) then  ! residuals
	write(*,'(a)') 
	write(*,'(a)') '<div class="ltgr">'
	write(*,'(a)') '# -------------------------------------------<br />'
	write(*,'(a)') '# Analysis of Residuals, Influence Statistics<br />'
	write(*,'(a)') '# -------------------------------------------<br />'
	write(*,'(a)') '</div>'
	txt1 = '<span class="b">fitted</span>(lm.out),<span class="b">resid</span>(lm.out)'
	write(*,'(a)') 'outr <- <span class="b">cbind</span>('//trim(txt1)//')<br />'
	txt1 = '<span class="b">rstudent</span>(lm.out),<span class="b">cooks.distance</span>(lm.out)'
	write(*,'(a)') 'outr <- <span class="b">cbind</span>(outr, '//trim(txt1)//')<br />'
	cm = 'outr <- <span class="b">cbind</span>(lm.out$model[c('
	cm = trim(cm) // '"<span class="r">'//trim(X(1))//'</span>", '
	if (nX > 1) then
		do i=2,nX
			cm = trim(cm) //  '"<span class="r">'//trim(X(i))//'</span>", '
		end do
	end if
	cm = trim(cm) // '"<span class="r">'//trim(Y)// '</span>")'
	write(*,'(a)') trim(cm)// '], outr)<br />'  
	write(*,'(a)') 'outr <- <span class="b">data.frame</span>(outr)<br />'
  write(*,'(a)') '<span class="b">names</span>(outr)['//c1//'] <- "fitted"<br />'
  write(*,'(a)') '<span class="b">names</span>(outr)['//c2//'] <- "residual"<br />'
  write(*,'(a)') '<span class="b">names</span>(outr)['//c3//'] <- "rstudent";<br />'
  write(*,'(a)') '<span class="b">names</span>(outr)['//c4//'] <- "cooks"<br />'
	write(*,'(a)') 'o <- <span class="b">order</span>(outr$cooks, decreasing=TRUE)<br />'
	write(*,'(a)') 'outr <- outr[o,]<br />'
	write(*,'(a)') 'outr[1:nrow(outr),]<br /><br />'
end if

if (pred_flag) then  ! confidence and prediction intervals
	write(*,'(a)') 
	write(*,'(a)') '<div class="ltgr">'
	write(*,'(a)') '# -------------------------------------------------------<br />'
	write(*,'(a)') '# Confidence and Prediction intervals from given X values<br />'
	write(*,'(a)') '# -------------------------------------------------------<br />'
	write(*,'(a)') '</div>'
	txt1 = 'ci <- <span class="b">data.frame</span>(<span class="b">predict</span>'
	write(*,'(a)') trim(txt1)//'(lm.out, interval="confidence"))<br />'
	txt1 = '<span class="b">data.frame</span>'
	txt2 = '(<span class="b">predict</span>(lm.out, interval="prediction"))'
	write(*,'(a)') 'pi <- suppressWarnings('//trim(txt1)//trim(txt2)//')<br />'
	cm = 'outr <- <span class="b">cbind</span>(lm.out$model[c('
	cm = trim(cm) // '"<span class="r">'//trim(X(1))//'</span>", '
	if (nX > 1) then
		do i=2,nX
			cm = trim(cm) //  '"<span class="r">'//trim(X(i))//'</span>", '
		end do
	end if
	cm = trim(cm) // '"<span class="r">'//trim(Y)// '</span>")'
	write(*,'(a)') trim(cm)// '],  ci,pi$lwr,pi$upr)<br />'  
	write(*,'(a)') 'outr <- <span class="b">data.frame</span>(outr)<br />'
	write(*,'(a)') 'o <- <span class="b">order</span>(outr[,'//c4//'], decreasing=FALSE)<br />'
	write(*,'(a)') 'outr <- outr[o,]<br />'
	write(*,'(a)') '<span class="b">names</span>(outr)['//c1//'] <- "fitted"<br />'
	write(*,'(a)') '<span class="b">names</span>(outr)['//c2//'] <- "ci:lwr"<br />'
	write(*,'(a)') '<span class="b">names</span>(outr)['//c3//'] <- "ci:upr"<br />'
	write(*,'(a)') '<span class="b">names</span>(outr)['//c4//'] <- "pi:lwr"<br />'
	write(*,'(a)') '<span class="b">names</span>(outr)['//c5//'] <- "pi:upr"<br />'
	write(*,'(a)') 'outr[1:nrow(outr),]<br />'
	write(*,'(a)') '<div class="ltgr">'
 	write(*,'(a)') '# Note: Predictions from current data, from which the model is estimated,<br />'
  write(*,'(a)') '#    refer to _future_ responses based on the collection of new data.<br />'
	write(*,'(a)') '</div>'
	write(*,'(a)') '<br />'
end if


if (best_flag) then ! best subset
	write(*,'(a)') '<div class="ltgr">'
	write(*,'(a)') '# ---------------------------------------------<br />'
	write(*,'(a)') '# R squared adjusted for all possible models<br />'
	write(*,'(a)') '# ---------------------------------------------<br />'
	write(*,'(a)') '# Warning:  This analysis is only of descriptive<br />'
	write(*,'(a)') '#   statistics, so does not literally generalize<br />'
	write(*,'(a)') '#   to the overall population.<br />'
	write(*,'(a)') '# ---------------------------------------------<br />'
	write(*,'(a)') '# one time only:  install.packages("leaps")<br />'
 	write(*,'(a)') '# leaps - find best subset based on adjusted R-squared<br />'
 	write(*,'(a)') '# create data frame X of just the predictor variables<br />'
	write(*,'(a)') '# modify output to get adjR2 column next to relevant model<br />'
	write(*,'(a)') '# sort by adjusted R2<br />'
	write(*,'(a)') '</div>'
	write(*,'(a)') '<span class="b">library</span>(<span class="b">leaps</span>)<br />'
	cm = 'X &lt;- <span class="b">data.frame</span>(<span class="b">cbind</span>(<span class="r">' 
	cm = trim(cm) // trim(X(1))// '</span>'
	if (nX > 1) then
		do i=2,nX
			cm = trim(cm) //',<span class="r">' // trim(X(i)) // '</span>'
		end do
	end if
	write(*,'(a)') trim(cm)// '))<br />'
	write(*,'(a)') 'leaps.out &lt;- <span class="b">leaps</span>(X, <span class="r">'//trim(Y)//'</span>, '
	write(*,'(a)') 'method="adjr2")<br />'
	txt1 = 'models &lt;- <span class="b">data.frame</span>'
	txt2 = '(leaps.out$which, row.names=NULL)<br />'
	write(*,'(a)') trim(txt1)//trim(txt2)
	write(*,'(a)') 'crit &lt;- <span class="b">data.frame</span>(leaps.out$adjr2)<br />'
	write(*,'(a)') '<span class="b">names</span>(crit) &lt;- "adjr2"<br />'
	write(*,'(a)') 'out &lt;- <span class="b">cbind</span>(models,crit)<br />'
	write(*,'(a)') 'sorted.out <- out[<span class="b">order</span>(out$adjr2, '
	write(*,'(a)') 'decreasing=TRUE),]<br />'
	write(*,'(a)') '<span class="b">names</span>(sorted.out) &lt;- c(<span class="b">names</span>(X),"Rsqadj")<br />'
	write(*,'(a)') 'sorted.out<br /><br />'
end if


if (nX == 1) then  ! plot scatterplot
	write(*,'(a)') 
	write(*,'(a)') '<div class="ltgr">'
	write(*,'(a)') '# -----------------------------------<br />'
	write(*,'(a)') '# Scatterplot and Regression Line<br />'
	write(*,'(a)') '# -----------------------------------<br />'
	if (pred_flag) then
		write(*,'(a)') '# Confidence and Prediction Intervals<br />'
		write(*,'(a)') '# -----------------------------------<br />'
	end if
	write(*,'(a)') '</div>'
	if (.NOT.pred_flag) then
		ctitle = '"Scatterplot and Regression Line"'
	else 
		ctitle = '"Regression Line, Confidence and Prediction Intervals"'
	end if
	txt1 = '<span class="b">plot</span>(<span class="r">'//trim(X(1))
	txt2 = '</span>, pch=19, cex=.8, col="gray70")<br />'
	write(*,'(a)') trim(txt1)//'</span>,<span class="r">'//trim(Y)//trim(txt2)
	txt1 = '<span class="b">title</span>'
	write(*,'(a)') trim(txt1)//'(main='//trim(ctitle)//')<br />'
	write(*,'(a)') '<span class="b">abline</span>(lm.out$coef)<br />'
	if (pred_flag) then
		txt1 = '<span class="b">lines</span>(<span class="r">'//trim(X(1))
		write(*,'(a)') trim(txt1)//'</span>, ci$lwr, col="lightsteelblue", lwd=2)<br />'
		write(*,'(a)') trim(txt1)//'</span>, ci$upr, col="lightsteelblue", lwd=2)<br />'
		write(*,'(a)') trim(txt1)//'</span>, pi$lwr, col="darksalmon", lwd=2)<br />'
		write(*,'(a)') trim(txt1)//'</span>, pi$upr, col="darksalmon", lwd=2)<br />'
	end if
else if (nX > 1) then  ! plot scatterplot matrix
	write(*,'(a)') '<div class="ltgr">'
	write(*,'(a)') '# ------------------<br />'
	write(*,'(a)') '# Scatterplot Matrix<br />'
	write(*,'(a)') '# ------------------<br />'
	write(*,'(a)') '</div>'
	txt1 = '<span class="b">pairs</span>(~ <span class="r">'//trim(Y)//'</span> + '
	write(*,'(a)') txt1//'<span class="r">'//trim(X(1))// '</span>'
	do i=2,nX
		cm = '+ <span class="r">' // trim(X(i)) // '</span>'
	end do
	write(*,'(a)') trim(cm)// ')<br /><br />'
end if

write(*,'(a)')  
write(*,'(a)') '</body>'
write(*,'(a)') '</html>'

end program regression
