## Who can help? [R for BE/BA]

Hi all,

some help needed.

Specifically, I am asking for help to make code work and not to make it smarter (achieving the same with less key strokes, using built-in alternatives etc) although there is generally a preferred focus on the latter here

Hopefully we will discuss

The idea is to construct the model from the bottoms up.

Another site for the objective function is: https://www.stt.msu.edu/users/pszhong/Lecture_16_Spring_2017.pdf (page 13+14), but this also did not quite get me to convergence.

So, this is what I have so far. I don't understand matrix algebra totally well, and I am at a loss of you say Hermitian, expectation-minimisation, Jacobian, etc.

What bothers me is that some years ago I wrote a working optimiser in C on basis of slides some some Danish guy called Waagepetersen, who walked the reader slowly through the process of optimising REML. I can't find my code now, it was on another computer, and I can't find his slides any more (there are some slides with his name out there, which I can find with google etc but they do not have that old recipe).

some help needed.

Specifically, I am asking for help to make code work and not to make it smarter (achieving the same with less key strokes, using built-in alternatives etc) although there is generally a preferred focus on the latter here

Hopefully we will discuss

**later**the need for additional or other variance components**once**we have made the optimiser work somehow (or**while**we are making it work, if such vc's are totally necessary). There are probably a ton of errors to iron out. There is at least one big one whose effects I a sufferingThe idea is to construct the model from the bottoms up.

- First I do a full-rank design matrix X for the fixed effects, without subjects

- Next I do the covariance matrix with just three variance components as discussed above.

- I generate some cheap starter guesses for the variance components

- I define an objective function for "a constant + log likelihood"

- I optimise the variance components, in REML the effect estimates will come along with them.

Another site for the objective function is: https://www.stt.msu.edu/users/pszhong/Lecture_16_Spring_2017.pdf (page 13+14), but this also did not quite get me to convergence.

So, this is what I have so far. I don't understand matrix algebra totally well, and I am at a loss of you say Hermitian, expectation-minimisation, Jacobian, etc.

What bothers me is that some years ago I wrote a working optimiser in C on basis of slides some some Danish guy called Waagepetersen, who walked the reader slowly through the process of optimising REML. I can't find my code now, it was on another computer, and I can't find his slides any more (there are some slides with his name out there, which I can find with google etc but they do not have that old recipe).

`rm(list=ls(all=TRUE))`

###### par vector is c(varT, varBR, varWR)

######

###### Section 1: Household things

######

CreateX=function(D)

{

##first the two treatments

TrtT=as.numeric(as.character(D$Trt)=="T")

TrtR=as.numeric(as.character(D$Trt)=="R")

X=cbind(TrtT, TrtR)

R= qr(X)$rank

## subjects: in a mixed model with subj as random

## we do not do subjects also as fixed, therefore they are #'ed away here

## for (s in unique(D$Subj))

## {

## v=as.numeric(D$Subj==s)

## #print(v)

## XX=data.frame(X, v)

## names(XX)[ncol(XX)] = paste("S", s, sep="")

## rnk=qr(XX)$rank

## if (rnk>R)

## {

## X=XX

## R=rnk

## }

## }

##now the Pers

for (p in unique(D$Per))

{

v=as.numeric(D$Per==p)

#print(v)

XX=data.frame(X, v)

names(XX)[ncol(XX)] = paste("P", p, sep="")

rnk=qr(XX)$rank

if (rnk>R)

{

X=XX

R=rnk

}

}

for (q in unique(D$Seq))

{

v=as.numeric(D$Seq==q)

#print(v)

XX=data.frame(X, v)

names(XX)[ncol(XX)] = paste("Q", q, sep="")

rnk=qr(XX)$rank

if (rnk>R)

{

X=XX

R=rnk

}

}

return(as.matrix(X))

}

Create.CovM=function(Params)

##block diagonal covariance matrix

{

varT=Params[1]

varBR=Params[2]

varWR=Params[3]

#varRT=Params[4]

#cat("Vars:", varT, varBR, varWR,"\n")

Nobs=length(D$Y)

V=matrix(0,ncol=Nobs, nrow=Nobs)

for (iRow in 1:Nobs)

for (iCol in 1:Nobs)

{

## the diagonal

if (iCol==iRow)

{

if (D$Trt[iRow]=="T") V[iRow,iCol]=V[iRow,iCol]+varT

if (D$Trt[iRow]=="R") V[iRow,iCol]=V[iRow,iCol]+varWR

}

## off diagonal

if (iCol!=iRow)

if (D$Subj[iRow]==D$Subj[iCol])

{

if (D$Trt[iCol]==D$Trt[iRow]) V[iRow,iCol]= V[iRow,iCol]+varBR

#if (D$Trt[iCol]!=D$Trt[iRow]) V[iRow,iCol]= V[iRow,iCol]+varRT

}

}

return(as.matrix(V))

}

######

###### Section 2: Matrix things

######

Obj.F12=function(Pars)

##this is line 3 of page 10 of:

##http://people.csail.mit.edu/xiuming/docs/tutorials/reml.pdf

{

CovM=Create.CovM(Pars)

A= -0.5*log(det(CovM))

B= -0.5*log(det(t(X) %*% solve(CovM) %*% X))

est.b = solve(t(X) %*% solve(CovM) %*% X) %*% t(X) %*% solve(CovM) %*% y

tmp= y - X %*% est.b

C=-0.5 *(t(tmp) %*% solve(CovM) %*% tmp)

return(A+B+C)

}

Some.Initial.Guesses=function(foo)

{

#guess variance within R

D1=subset(D, D$Trt=="R")

m=lm(log(Y)~factor(Subj)+factor(Seq)+factor(Per) , data=D1)

#print(anova(m))

varWR= anova(m)["Residuals", "Mean Sq"]

m=lm(log(Y)~factor(Subj)+factor(Seq)+factor(Per)+factor(Trt) , data=D)

varBR= anova(m)["factor(Subj)", "Mean Sq"]

##total var T, cheaply

##actually, this may be much better: varT=varBR+varWR

D2=subset(D, D$Trt=="T")

varT=var(log(D2$Y))

varT=varBR+varWR

L=list(varT=varT, varBR=varBR, varWR=varWR)

return(L)

}

####

#### section 3: execution

####

D=read.csv("EMAII.a.csv" , header=T, sep="\t") ##public = easier

X=CreateX(D) ##public = easier

y=log(D$Y) ##public = easier

Check.Surf=function(i, x1,x2, N)

## Checks the log likelihood surface of the i'th variance component

## (given the other variance components) for the interval var = x1 to x2

## with a granularity of N points.

{

L=Some.Initial.Guesses(0)

parIni=c(L$varT, L$varBR, L$varWR )

x=x1;

vX=NULL

vY=NULL

dx=(x2-x1)/N

for (j in 1:N)

{

vX=c(vX, x)

parIni[i]=x

vY=c(vY, Obj.F12(parIni))

x=x+dx

}

plot(vX, vY)

}

MyREML=function(foo.bar)

## ideally this would work, but it doesn't right now

## also, I have a feeling the surface may be quite rugged,

## so I may switch to a homecooked optimiser

## once I have ironed out the singularity issue, see below.

{

L=Some.Initial.Guesses(0)

print(L)

parIni=c(L$varT, L$varBR, L$varWR)

F=optim(par=parIni, fn=Obj.F12, #method="L-BFGS-B",lower=c(0,0,0),

control=list(reltol=1e-10, trace=T))

print(F)

}

# if ti worked it'd be called through

# MyREML(0)

# however:

Check.Surf(1, 0.01, 0.2, 50)

## ouch, that's a singularity.

## so something is wrong in the covariance matrix

## or in the "const + likelihood function"

## or elsewhere. Any ideas?

#

#

#

—

I could be wrong, but...

Best regards,

ElMaestro

R's base package has 274 reserved words and operators, along with 1761 functions. I can use 18 of them (about 14 of them properly). I believe this makes me the Donald Trump of programming.

I could be wrong, but...

Best regards,

ElMaestro

R's base package has 274 reserved words and operators, along with 1761 functions. I can use 18 of them (about 14 of them properly). I believe this makes me the Donald Trump of programming.

### Complete thread:

- Semireplicated + REML in R ElMaestro 2020-07-10 13:19 [R for BE/BA]
- Semireplicated + REML in R Helmut 2020-07-10 19:03
- Semireplicated + REML in R ElMaestro 2020-07-10 19:24
- Avoid partial replicate designs, pleeeze! Helmut 2020-07-11 11:57
- Avoid partial replicate designs, pleeeze! ElMaestro 2020-07-11 14:43
- Braveheart! Helmut 2020-07-11 15:38
- Who can help?ElMaestro 2020-07-12 12:28
- Who can help? ElMaestro 2020-07-12 13:26
- Update II ElMaestro 2020-07-12 21:36
- Update III ElMaestro 2020-07-12 21:46

- Final update today ElMaestro 2020-07-12 22:27
- Medium rare. Helmut 2020-07-13 13:52
- took just 52 hrs to do it :-) ElMaestro 2020-07-13 14:18
- Will take much more hours still… Helmut 2020-07-13 15:34

- Negative determinant ElMaestro 2020-07-14 03:22
- Are we loosers? Helmut 2020-07-14 13:58
- "we"? Loosers? ElMaestro 2020-07-14 15:07
- Misunderstanding? Helmut 2020-07-14 15:32

- "we"? Loosers? ElMaestro 2020-07-14 15:07

- Are we loosers? Helmut 2020-07-14 13:58

- took just 52 hrs to do it :-) ElMaestro 2020-07-13 14:18

- Medium rare. Helmut 2020-07-13 13:52

- Braveheart! ElMaestro 2020-07-13 10:13
- Braveheart! Helmut 2020-07-13 14:16
- Braveheart! PharmCat 2020-07-15 14:19
- Braveheart! ElMaestro 2020-08-02 17:39

- Who can help?ElMaestro 2020-07-12 12:28

- Braveheart! Helmut 2020-07-11 15:38

- Avoid partial replicate designs, pleeeze! ElMaestro 2020-07-11 14:43

- Avoid partial replicate designs, pleeeze! Helmut 2020-07-11 11:57

- Semireplicated + REML in R ElMaestro 2020-07-10 19:24
- We were all blind (except Detlew) Helmut 2020-07-15 14:27
- It is the opposite way around for me ElMaestro 2020-07-15 16:25
- Desultory thoughts Helmut 2020-07-15 17:33
- FDA RSABE is ISC d_labes 2020-07-15 18:13
- FDA RSABE is ISC Helmut 2020-07-16 11:11

- Desultory thoughts ElMaestro 2020-07-15 23:06
- Desultory thoughts Helmut 2020-07-16 10:59

- FDA RSABE is ISC d_labes 2020-07-15 18:13

- Desultory thoughts Helmut 2020-07-15 17:33
- Phoenix - which template? mittyri 2020-07-19 00:42
- FDA RSABE Project template_ v1.4.phxproj Helmut 2020-07-19 01:45

- It is the opposite way around for me ElMaestro 2020-07-15 16:25
- "By popular demand": likelihood ElMaestro 2020-07-24 10:07
- And by the way.... ElMaestro 2020-07-24 12:52
- And by the way.... PharmCat 2020-08-03 14:24
- Not understood ElMaestro 2020-08-03 22:55
- Not understood PharmCat 2020-08-05 01:41
- Not understood ElMaestro 2020-08-05 08:13
- Not understood PharmCat 2020-08-05 16:37
- Open issues ElMaestro 2020-08-06 21:11
- Open issues PharmCat 2020-08-07 00:02
- Open issues ElMaestro 2020-08-07 07:49
- Open issues PharmCat 2020-08-07 11:29
- Still can't make it work ElMaestro 2020-08-07 13:08
- Still can't make it work PharmCat 2020-08-07 15:42
- Still can't make it work ElMaestro 2020-08-07 16:44
- Still can't make it work PharmCat 2020-08-07 18:14
- Still can't make it work ElMaestro 2020-08-07 18:23
- And now it works ElMaestro 2020-08-07 21:31
- Still can't make it work PharmCat 2020-08-08 01:08
- Speed improvement ElMaestro 2020-08-08 12:25
- Speed improvement PharmCat 2020-08-08 17:27
- Speed improvement ElMaestro 2020-08-08 18:10
- Speed improvement PharmCat 2020-08-09 18:22
- Some tests... PharmCat 2020-08-10 11:48

- Speed improvement ElMaestro 2020-08-08 12:25

- Still can't make it work ElMaestro 2020-08-07 18:23

- Still can't make it work PharmCat 2020-08-07 18:14

- Still can't make it work ElMaestro 2020-08-07 16:44

- Still can't make it work PharmCat 2020-08-07 15:42

- Still can't make it work ElMaestro 2020-08-07 13:08

- Open issues PharmCat 2020-08-07 11:29

- Open issues ElMaestro 2020-08-07 07:49

- Open issues PharmCat 2020-08-07 00:02

- Open issues ElMaestro 2020-08-06 21:11

- Not understood PharmCat 2020-08-05 16:37

- Not understood ElMaestro 2020-08-05 08:13

- Not understood PharmCat 2020-08-05 01:41

- Not understood ElMaestro 2020-08-03 22:55

- And by the way.... PharmCat 2020-08-03 14:24

- And by the way.... ElMaestro 2020-07-24 12:52

- Semireplicated + REML in R Helmut 2020-07-10 19:03