Re: [R] ipf function in R

2008-03-06 Thread GREGOR Brian J
Chandra,

While I can't advise on how to use the ipf function in the cat package,
I can offer the following function that we use here to balance rebalance
arrays with new marginal totals.

#ipf.R
#Function to iteratively proportionally fit a multidimensional array
#IPF also known as Fratar method, Furness method, raking and two/three
dimensional balancing.
#This method of matrix balancing is multiplicative since the margin
factors (coefficients)
#are multiplied by the seed array to yield the balanced array. 
#Ben Stabler, [EMAIL PROTECTED], 9.30.2003
#Brian Gregor, [EMAIL PROTECTED], 2002

#inputs:
#1) margins_ - a list of margin values with each component equal to a
margin
#Example: row margins of 210 and 300, column margins of 140 and 370
#and third dimension margins of 170 and 340.  
#[[1]] [,1] [,2] [,3]
#210, 300
#[[2]]
#140, 370
#[[3]]
#170, 340
#2) seedAry - a multi-dimensional array used as the seed for the IPF
#3) iteration counter (default to 100)
#4) closure criteria (default to 0.001)

#For more info on IPF see:
#Beckmann, R., Baggerly, K. and McKay, M. (1996). Creating Synthetic
Baseline Populations.
#   Transportation Research 30A(6), 415-435.
#Inro. (1996). Algorithms. EMME/2 User's Manual. Section 6.
#~~~
~~~

ipf - function(Margins_, seedAry, maxiter=100, closure=0.001) {
#Check to see if the sum of each margin is equal
MarginSums. - unlist(lapply(Margins_, sum))
if(any(MarginSums. != MarginSums.[1])) warning(sum of each margin
not equal)

#Replace margin values of zero with 0.001
Margins_ - lapply(Margins_, function(x) {
  if(any(x == 0)) warning(zeros in marginsMtx replaced with
0.001) 
  x[x == 0] - 0.001
  x
  })

#Check to see if number of dimensions in seed array equals the
number of
#margins specified in the marginsMtx
numMargins - length(dim(seedAry))
if(length(Margins_) != numMargins) {
stop(number of margins in marginsMtx not equal to number of
margins in seedAry)
}

#Set initial values
resultAry - seedAry
iter - 0
marginChecks - rep(1, numMargins)
margins - seq(1, numMargins)

#Iteratively proportion margins until closure or iteration criteria
are met
while((any(marginChecks  closure))  (iter  maxiter)) {
for(margin in margins) {
marginTotal - apply(resultAry, margin, sum)
marginCoeff - Margins_[[margin]]/marginTotal
marginCoeff[is.infinite(marginCoeff)] - 0
resultAry - sweep(resultAry, margin, marginCoeff, *)
marginChecks[margin] - sum(abs(1 - marginCoeff))
}
iter - iter + 1
}

#If IPF stopped due to number of iterations then output info
if(iter == maxiter) cat(IPF stopped due to number of iterations\n)

#Return balanced array
resultAry
}


Brian Gregor, P.E.
Transportation Planning Analysis Unit
Oregon Department of Transportation
[EMAIL PROTECTED]
(503) 986-4120

 
 Message: 143
 Date: Wed, 05 Mar 2008 18:14:28 +1100
 From: Chandra Shah [EMAIL PROTECTED]
 Subject: [R] ipf function in R
 To: r-help@r-project.org
 Message-ID: [EMAIL PROTECTED]
 Content-Type: text/plain; charset=ISO-8859-1; format=flowed
 
 Hi
 I have a 3 x 2 contingency table:
 10 20
 30 40
 50 60
 I want to update the frequencies to new marginal totals:
 100 130
 40 80 110
 I want to use  the ipf (iterative proportional fitting) 
 function which 
 is apparently in the cat package.
 Can somebody please advice me how to input this data and 
 invoke ipf in R 
 to obtain an updated contingency table?
 Thanks.
 By the way I am quite new to R.
 
 -- 
  
 
 Dr Chandra Shah
 Senior Research Fellow
 Monash University-ACER Centre for the Economics of Education 
 and Training
 Faculty of Education, Building 6,
 Monash University
 Victoria
 Australia 3800
 Tel. +61 3 9905 2787
 Fax +61 3 9905 9184
 
 www.education.monash.edu.au/centres/ceet
 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] ipf function in R

2008-03-05 Thread Ted Harding
On 05-Mar-08 07:14:28, Chandra Shah wrote:
 Hi
 I have a 3 x 2 contingency table:
 10 20
 30 40
 50 60
 I want to update the frequencies to new marginal totals:
 100 130
 40 80 110
 I want to use  the ipf (iterative proportional fitting) function
 which is apparently in the cat package.
 Can somebody please advice me how to input this data and invoke
 ipf in R to obtain an updated contingency table?
 Thanks.
 By the way I am quite new to R.

It is not clear from your description what you really want to do.
Hence, not clear whether you should be using ipf() to do it!

The purpose of ipf() in the 'cat' package is very simple:
You suppply a complete contingency table, and specify the
order of interactions you want to include in the fit (using
the margins parameter in ipf()); and then the value of ipf()
is a table of the fitted (expected) values corresponding to
a model of the type you have specified, fitted by maximum
likelihood.

Your task seems to be different. If I understand your statement.
it looks as though you want to find a,b,c,d such that

  10+a  20+c | 40
  30+b  40+d | 80
  5060   |110
-+---
 100   130   |230

or, since there are row and column constraints on a,b,c,d:

  10+a  30-a | 40
  40-a  40+a | 80
  5060   |110
-+---
 100   130   |230

so then the task would be to find a.

But now it depends on what you would mean by find a.

You could treat 'a' as an unknown parameter, and fit it
by maximum likelihood under various assumptions about the
dependence between rows and columns (e.g. that the log
odds-ratios have particular values).

Or perhaps (and here is where ipf() may come in) for a
given choice of 'a' you could submit the resulting table
to ipf() and compute chisq from the resulting fitted
values and the given observed values; then seek the
value of 'a' which minimises chisq.

However, it's not really possible to advise in more detail
without being sure of what you are trying to do. Above,
I am only guessing!

Best wishes,
Ted.


E-Mail: (Ted Harding) [EMAIL PROTECTED]
Fax-to-email: +44 (0)870 094 0861
Date: 05-Mar-08   Time: 10:26:57
-- XFMail --

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] ipf function in R

2008-03-05 Thread Charles C. Berry
On Wed, 5 Mar 2008, Chandra Shah wrote:

 Hi
 I have a 3 x 2 contingency table:
 10 20
 30 40
 50 60
 I want to update the frequencies to new marginal totals:
 100 130
 40 80 110
 I want to use  the ipf (iterative proportional fitting) function which
 is apparently in the cat package.
 Can somebody please advice me how to input this data and invoke ipf in R
 to obtain an updated contingency table?

I'd use loglin()

newtab -
loglin( rowmarg%o%colmarg/sum(colmarg),
margin=list(1,2),
start=tab, fit=TRUE )$fit


with rowmarg and colmarg set to your updated marginals.

As for inputting the data, if this is all you have, type it in at the 
command line. See

?matrix
?c

and note this:

 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


HTH,

Chuck


 Thanks.
 By the way I am quite new to R.

 -- 


 Dr Chandra Shah
 Senior Research Fellow
 Monash University-ACER Centre for the Economics of Education and Training
 Faculty of Education, Building 6,
 Monash University
 Victoria
 Australia 3800
 Tel. +61 3 9905 2787
 Fax +61 3 9905 9184

 www.education.monash.edu.au/centres/ceet

 __
 R-help@r-project.org mailing list
 https://stat.ethz.ch/mailman/listinfo/r-help
 PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
 and provide commented, minimal, self-contained, reproducible code.


Charles C. Berry(858) 534-2098
 Dept of Family/Preventive Medicine
E mailto:[EMAIL PROTECTED]  UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.