runpool<-function(prefix,age="teen",data=workdat1,dat3){
  #Create a list of all of the variables for which we want listwise
  #deletion, (since some of them won't be in this particular
  #regression, the next few steps deal with doing the listwise
  #deletion on the relevant part of the pooled dataset --- in an
  #extraordinarily clunky manner --- taking the data apart into two
  #pieces, doing listwise deletion on each one, putting them back
  #together.Yuck!) ANYWAY, the next step chooses those variables in
  #the dataset whose names contain the prefix given above.
  
  allvars<-c("y65id","pan",names(data)[grep(c(prefix),names(data))])
  #Is this regression for the teens or the twenty-year-olds?
  agen<-switch(age,"teen"=1,"twen"=2)
  #Select the variables using the regression.
  ycol<-allvars[grep(c(paste("kid",agen,sep="",coll="")),allvars)]
  xcol<-allvars[grep(c(paste("par",agen,sep="",coll="")),allvars)]
  datcols<-c("y65id","pan",ycol,xcol)
  dat1<-data[data$pan==0,datcols]
  dat2<-na.omit(data[data$pan==1,allvars])
  dat3<-rbind(dat1,dat2[,datcols])
#  attach(dat3)
  #Create the regression formula as a text object, and parse it.
  form1<-parse(text = paste(ycol,"~",xcol,"*pan",sep = "", collapse = ""))
  lm1<-lm(formula=eval(form1),data=dat3,na.action=na.omit,model=T)
  summary(lm1)
  N1<-length(dat3$pan==0)
  N2<-length(dat3$pan==1)
  #Finite sample correction term, the 2
  #in the denom refers to the fact that, in separate regressions, 2
  #coeffcients are estimated
  fincorr<-function(n) (n-1)/(n-2)
  #Create the weight using the variances of the residuals from the previous regression                  
  wt1<-ifelse(dat3$pan==0,  
              var(lm1$resid[lm1$model$pan==0])*fincorr(N1),
              var(lm1$resid[lm1$model$pan==1])*fincorr(N2))
  #Run WLS
  lm2<-update(lm1,weight=1/wt1)
#  detach(dat3)
  summary(lm2)
}

reportpool<-function(x){
  #Calculates relevant effects and takes the p-value for the
  #interaction term off of the summary object created above. Writes
  #out a 1x3 vector.
  junkie<-c(x$coef[2,1]+x$coef[4,1],x$coef[2,1],x$coef[4,4])
  names(junkie)<-c("G1G2","G2G3","Pval")
  junkie}

waldtest<-function(olsobj){
#Wald Statistic after OLS (I haven't tested it after glm)
  vcov<-summary(olsobj)$sigma^2*summary(olsobj)$cov.unscaled
  #const<-t(c(0,0,1,1)) #tests for b3+b4=0
  #q<-0 
  #(see Green 337 for hints about setting up different constraint matrices)
  #const<-matrix(c(0,0,1,0,0,0,0,1),2,4,byrow=T) #tests for b3=0 &b4=0
                                        #(can also do anova(lm1,lm2,test="Chisq") on two models)
  #q<-c(0,0)
  beta<-as.matrix(olsobj$coef)
  wald<-t(const%*%beta-q)%*%solve(const%*%vcov%*%t(const))%*%(const%*%beta-q)
  pval<-1-pchisq(wald,nrow(const))
  outjunk<-c(wald,nrow(const),pval)
  names(outjunk)<-c("Wald","df","P Chisq")
  outjunk
}