data skin;
    retain id 0;
    infile 'skin.txt' firstobs=2 expandtabs;
    input cases towncd age $ population;
    town="Minneapolis";
    if  towncd=1 then town="Dallas";
    logpop=log(population);
	rate=cases/population;
    id=id+1;

*%let ginv=sweep;
%include "glmm.sas";

proc iml;
    load _all_;
    use skin;
    read all;

    nobs=nrow(cases);
    x=j(nobs,1,1)||design(town)||design(age);
    
    effname={"Intercept"}||(unique(town))||(unique(age));
    effname=effname`;
    print x[colname=effname] ;
    start link(y,eta,hrh,hry,wt,log_like,mean,r);
        BIG=20;
        eta=BIG*(eta>BIG)-BIG*(eta<-BIG)+eta#(abs(eta)<=BIG);
       call poi_link(y,eta,hrh,hry,wt,log_like,mean,r);
    finish link;

    eta=j(nobs,1,-5);
    call glm(x,cases,eta,beta,lhs,rhs,log(population),log_like,iter,delt,200,1.e-8,0);
    cov=sweep(lhs);
    stderr=sqrt(vecdiag(sweep(lhs)));
    print "Final";
    print effname beta stderr log_like delt iter ;
    dfe=nobs-round(trace(ginv(lhs)*lhs));
    etay=log(cases/population);
    offset=log(population);
    call link(cases,eta,hrh,hry,offset,log_like,mean,r);
    call link(cases,etay,hrh,hry,offset,dev,meany,ry);
    deviance=2*(dev-log_like);
    devscale=sqrt(deviance/dfe);
    pearscale=sqrt((cases-mean)`*inv(r)*(cases-mean)/dfe);
    print devscale pearscale dfe ;
        
    K=({0 1 -1 }`)//j(ncol(unique(age)),1,0);
    
    label="City";
    call wald(lhs,beta,k,TS,df,pvalue,label);

    eta=j(nobs,1,2.5);
    call glm(x,cases,eta,beta,lhs,rhs,log(population),log_liker,iter,delt,200,1.e-8,k);
    stderr=sqrt(vecdiag(sweep(lhs)));
    print "Final";
    print effname beta stderr log_liker delt iter ;
    dfer=nobs-round(trace(ginv(lhs)*lhs));
    
    cov=k`*sweep(lhs)*k;
    df=dfer-dfe;
    G=2*(log_like-log_liker);
    pval=1-probchi(G,df);
    print label G df pval;
    
quit;
    
proc genmod data=skin;
    class town age;
    model cases =town age/dist=poi type3 type1 offset=logpop;
    estimate 'City M-D' town -1 1/exp;
    run;
proc genmod data=skin;
    class town age;
    model cases=town age/dist=poi type3 type1 dscale offset=logpop;
    estimate 'City M-D' town -1 1/exp;
    run;
proc genmod data=skin;
    class town age;
	weight population;
    model rate=town age/dist=poi type3 type1 dscale;
    estimate 'City M-D' town -1 1/exp;
    run;
proc glmmod data=skin outdesign=skin_design;
    class town age;
    model cases population id=town age;
data skin;
    merge skin skin_design;
proc sort data=skin;
    by id;
proc nlmixed data=skin;
    parms int -5 C1-c2=0 age1-age8=0 sigu 0;
    array eff{11} int  C1-C2  age1-age8 ;
    array X{11} col1-col11;
    eta=u;
    do i=1 to 11;
        eta=eta+eff{i}*X{i};
        end;
    mu=exp(eta)*population;
    model cases~poisson(mu);
    random u~normal(0,exp(sigu)) subject=id;
    estimate 'sig2u' exp(sigu);
    estimate 'Perc Red' 100*(1-exp(c2-c1));
run;