Inverse Probability Censoring Weighting Estimation

https://documentation.sas.com/?docsetId=statug&docsetTarget=statug_rmstreg_details13.htm&docsetVersion=15.1&locale=en


https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5379776/

BLUP for mixed model

http://rpubs.com/enwuliu/548305

This is R code for the example from the paper "That Blup Is a Good Thing: The Estimation of Random Effects.".

Robinson, George K. "That Blup Is a Good Thing: The Estimation of Random Effects." Statistical Science 6, no. 1 (1991): 15-32.

In the book Generalized Additive Models An Introduction with R, the second version, page 81 the function llm that solved the linear mixed model also used Henderson's equations but did not mention the equations.





mixed_model_blup.R
enliu
2019-11-10
Y<-matrix(c(110,100,110,100,100,110,110,100,100),nrow=9,ncol=1)
X<-matrix(c(1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,1),nrow=9,ncol=3)

Z<-matrix(c(1,0,0,0,0,0,0,0,0,
            0,0,1,0,0,0,0,0,0,
            0,0,0,0,0,1,1,0,0,
            0,1,0,1,1,0,0,1,1),nrow=9,ncol=4)
R<-diag(1,nrow=9,ncol=9)
G<-diag(0.1,nrow=4,ncol=4)
A<-t(X)%*%solve(R)%*%X
B<-t(X)%*%solve(R)%*%Z
C<-t(Z)%*%solve(R)%*%X
D<-t(Z)%*%solve(R)%*%Z+solve(G)

E<-cbind(A,B)
F<-cbind(C,D)

M<-rbind(E,F)

t(X)%*%solve(R)%*%Y
##      [,1]
## [1,]  210
## [2,]  310
## [3,]  420
t(Z)%*%solve(R)%*%Y
##      [,1]
## [1,]  110
## [2,]  110
## [3,]  220
## [4,]  500
N<-rbind(t(X)%*%solve(R)%*%Y,t(Z)%*%solve(R)%*%Y)

solve(M,N)
##             [,1]
## [1,] 105.6386574
## [2,] 104.2757378
## [3,] 105.4584366
## [4,]   0.3964857
## [5,]   0.5203875
## [6,]   0.7569272
## [7,]  -1.6738004



/*********Another method from the Generalized Additive Models An Introduction with R, the second version, p81


llm <- function(theta,X,Z,y) {
  ## untransform parameters...
  sigma.b <- exp(theta[1])
  sigma <- exp(theta[2])
  ## extract dimensions...
  n <- length(y); pr <- ncol(Z); pf <- ncol(X)
  ## obtain \hat \beta, \hat b...
  X1 <- cbind(X,Z)
  ipsi <- c(rep(0,pf),rep(1/sigma.b^2,pr))
 
  b1 <- solve(crossprod(X1)/sigma^2+diag(ipsi),t(X1)%*%y/sigma^2)
 
  ## compute log|Z'Z/sigma^2 + I/sigma.b^2|...
 
  ldet <- sum(log(diag(chol(crossprod(Z)/sigma^2 +
                              diag(ipsi[-(1:pf)])))))
  ## compute log profile likelihood...
  l <- (-sum((y-X1%*%b1)^2)/sigma^2 - sum(b1^2*ipsi) -n*log(sigma^2) - pr*log(sigma.b^2) - 2*ldet - n*log(2*pi))/2
  attr(l,"b") <- as.numeric(b1) ## return \hat beta and \hat b
  -l
}

y<-matrix(c(110,100,110,100,100,110,110,100,100),nrow=9,ncol=1)
x<-matrix(c(1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,1),nrow=9,ncol=3)

z<-matrix(c(1,0,0,0,0,0,0,0,0,
            0,0,1,0,0,0,0,0,0,
            0,0,0,0,0,1,1,0,0,
            0,1,0,1,1,0,0,1,1),nrow=9,ncol=4)


X1 <- cbind(x,z)

crossprod(X1)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    2    0    0    1    0    0    1
## [2,]    0    3    0    0    1    0    2
## [3,]    0    0    4    0    0    2    2
## [4,]    1    0    0    1    0    0    0
## [5,]    0    1    0    0    1    0    0
## [6,]    0    0    2    0    0    2    0
## [7,]    1    2    2    0    0    0    5
n <- length(y); pr <- ncol(z); pf <- ncol(x)

theta<-c(log(sqrt(0.1)),log(1))

sigma.b <- exp(theta[1])
sigma <- exp(theta[2])

ipsi <- c(rep(0,pf),rep(1/sigma.b^2,pr))


llm(theta,x,z,y)
## [1] 92.34967
## attr(,"b")
## [1] 105.6386574 104.2757378 105.4584366   0.3964857   0.5203875   0.7569272
## [7]  -1.6738004

crossprod function

http://www.endmemo.com/program/R/crossprod.php

R crossprod Function

crossprod() function returns matrix cross-product.

crossprod(x, y = NULL)
tcrossprod(x, y = NULL)

x: numeric matrix
y: numeric matrix, if y=NULL, y is the same as x
...

> x <- matrix(1:9,3,3)
> x

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> crossprod(x)

     [,1] [,2] [,3]
[1,]   14   32   50
[2,]   32   77  122
[3,]   50  122  194

> tcrossprod(x)

     [,1] [,2] [,3]
[1,]   66   78   90
[2,]   78   93  108
[3,]   90  108  126

R Tutorials

R Data Types

Loop, Condition Statements

Plotting and Graphics

String Manipulations

Math Functions

Matrix Manipulations

Read & Write Data

Statistical Analysis

All Functions List

:: Popular ::

»	R PCH Symbols
»	R Color Names
»	R Regular Expression
»	R tapply Function
»	R String Functions
»	R Plot Function
»	R Builtin Datasets List
	Python Tutorials
	HTML Tutorials
	JavaScript Tutorials
	Statistics
	News, Events Worldwide
	Unit Conversions
	Top Visited Websites Directory

Schur complement

https://en.wikipedia.org/wiki/Schur_complement

Output variance -covariance matrix from Genmod

proc genmod data=outmi;
   model Oxygen= RunTime RunPulse/covb;
   ods output ParameterEstimates=gmparms
              ParmInfo=gmpinfo
              CovB=gmcovb;
run;

Mean and Variance of epxonential family

https://math.stackexchange.com/a/1634096/166263

Piecewise linear basis, tent function, B-spline in GAM

require(gamair)
data("engine")
engine
attach(engine)
plot(size,wear,xlab='Engine capacity',ylab='Wear index')


tf<-function(x,xj,j){
  dj<-xj*0;
  dj[j]<-1
  approx(xj,dj,x)$y # approx return two subject, only return y here

}
tf.X<-function(x,xj){
  ##tent function basis matrix given data x
  ## and knot sequnce xj
  nk<-length(xj);
  n<-length(x)
  X<-matrix(NA,n,nk)
  for(j in 1:nk){
    X[,j]<-tf(x,xj,j)
  }
  X
}
sj<-seq(min(size),max(size),length=6)
sj
X<-tf.X(size,sj)
X

b<-lm(wear~X-1)
b
s<-seq(min(size),max(size),length=200)
Xp<-tf.X(s,sj)
plot(size,wear)
lines(s,Xp%*%coef(b))

Xp
library(splines)
B=bs(size,knots=sj,Boundary.knots=c(1.420,2.980),degree=1)
B

b2<-lm(wear~B-1)
b2

s2<-seq(min(size),max(size),length=200)
s2
Xp2<-bs(s2,knots=sj,Boundary.knots=c(1.420,2.980),degre=1)
Xp2[,1:6]

plot(size,wear)
lines(s2,Xp2[,1:6]%*%coef(b))

B=bs(s,knots=sj,Boundary.knots=c(1.420,2.980),degree =1)

matplot(s,B,type="l")