2020年10月24日星期六

Iterative re-weighted least squares(IRLS) for Poisson regression with identity link and canonical log link

 

y<-c(2,3,6,7,8,9,10,12,15)
x<-matrix(c(1,1,1,1,1,1,1,1,1,-1, -1, 0, 0, 0, 0, 1, 1, 1), nrow=9,ncol=2)
b.init<-c(7,5)
x1<-c(-1, -1, 0, 0, 0, 0, 1, 1, 1)
data_4.3<-data.frame(y,x1)

poisson_IRLS<-function(x,y,b.init,tol=1e-8){
  change <- Inf
  b.old <- b.init
  while(change > tol) {
    eta <- x %*% b.old  # linear predictor
    w<-diag(1/as.vector(eta),nrow=nrow(x),ncol=nrow(x))
    z<-y
    b.new<-solve(t(x)%*%w%*%x)%*%t(x)%*%w%*%z
    change <- sqrt(sum((b.new - b.old)^2))
    b.old <- b.new
  }
  b.new
}

poisson_IRLS(x,y,b.init)
##          [,1]
## [1,] 7.451633
## [2,] 4.935300
m1<-glm(y ~ x1, family=poisson(link="identity"),data=data_4.3)
summary(m1)
## 
## Call:
## glm(formula = y ~ x1, family = poisson(link = "identity"), data = data_4.3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7019  -0.3377  -0.1105   0.2958   0.7184  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.4516     0.8841   8.428  < 2e-16 ***
## x1            4.9353     1.0892   4.531 5.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.4206  on 8  degrees of freedom
## Residual deviance:  1.8947  on 7  degrees of freedom
## AIC: 40.008
## 
## Number of Fisher Scoring iterations: 3
inv_link<-function(x){exp(x)}

poisson_IRLS_log_link<-function(x,y,b.init,inv_link,tol=1e-8){
  change <- Inf
  b.old <- b.init
  while(change > tol){
    eta <- x %*% b.old  # linear predictor
    y.hat<-inv_link(eta)  #y.hat is mu in fact
    w<-diag(as.vector(y.hat),nrow=nrow(x),ncol=nrow(x))
    z<-solve(w)%*%(y-y.hat)
    b.new<-b.old+solve(t(x)%*%w%*%x)%*%t(x)%*%w%*%z
    change <- sqrt(sum((b.new - b.old)^2))
    b.old <- b.new
  }
  b.new
}

poisson_IRLS_log_link(x,y,rep(1,2),inv_link)
##           [,1]
## [1,] 1.8892720
## [2,] 0.6697856
m2<-glm(y ~ x1, family=poisson(link="identity"),data=data_4.3)
summary(m2)
## 
## Call:
## glm(formula = y ~ x1, family = poisson(link = "identity"), data = data_4.3)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.7019  -0.3377  -0.1105   0.2958   0.7184  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.4516     0.8841   8.428  < 2e-16 ***
## x1            4.9353     1.0892   4.531 5.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.4206  on 8  degrees of freedom
## Residual deviance:  1.8947  on 7  degrees of freedom
## AIC: 40.008
## 
## Number of Fisher Scoring iterations: 3

2020年6月6日星期六

SAS Macro functions – %eval() and %sysevalf()

The %EVAL function evaluates integer arithmetic or logical expressions. %EVAL operates by converting its argument from a character value to a numeric or logical expression. Then, it performs the evaluation. Finally, %EVAL converts the result back to a character value and returns that value.
Operands that contain a period character cause an error when they are part of an integer arithmetic expression. The following examples show correct and incorrect usage, respectively:
%let d=%eval(10+20);     /* Correct usage */
%let d=%eval(10.0+20.0); /* Incorrect usage */
When %EVAL encounters a value containing a period, it displays an error message about finding a character operand where a numeric operand is required.
%let a=1+2;
%let b=10*3;
%let c=5/3;
%let eval_a=%eval(&a);
%let eval_b=%eval(&b);
%let eval_c=%eval(&c);
 
%put &a is &eval_a;
%put &b is &eval_b;
%put &c is &eval_c;

The %SYSEVALF function evaluates arithmetic and logical expressions using floating-point arithmetic and returns a value that is formatted using the BEST32. format. The result of the evaluation is always text. % SYSEVALF is the only macro function that can evaluate logical expressions that contain floating point or missing values. Specifying a conversion type can prevent problems when %SYSEVALF returns missing or floating point values to macro expressions or macro variables that are used in other macro expressions that require an integer value.
Syntax: %SYSEVALF(expression<, conversion-type>)
Taking an example…
%let a=100;
%let b=101.590;
%let x=%sysevalf(&a+&b);
%let y=%sysevalf(&a +&b, boolean);
%let z=%sysevalf(&a +&b, ceil);
%let v=%sysevalf(&a +&b, floor);
%let w=%sysevalf(&a +&b, int);
 
%put The result with SYSEVALF is: &x;
%put  The BOOLEAN value is: &y;
%put  The CEIL value is: &z;
%put  The FLOOR value is: &v;
%put  The INTEGER value is: &w;

http://technico.qnownow.com/sas-macro-functions-eval-and-sysevalf-2/

2020年6月2日星期二

parmbuff (or Pbuff) in SAS macro



%macro demo(a=1, b=2)/parmbuff;
%put buffer holds |&syspbuff|;
%put &=a;
%put &=b;
%mend demo;
%demo(a=aa)

%demo(a=silly, d=unknown)