set.seed(1)
data <- data.frame(Used = rep(c(1,0,0,0),1250),
               Open = round(runif(5000,0,50),0),
               Activity = rep(sample(runif(24,.5,1.75),1250, replace=T), each=4),
               Strata = rep(1:1250,each=4))

mod <- clogit(Used ~ Open + I(Open*Activity) + strata(Strata),data=data)

newdata <- data.frame(Open = seq(0,50,1),
                  Activity = rep(max(data$Activity),51))

fit<-predict(mod,newdata=newdata,type = "expected")

coef<-data.frame(coef = summary(mod)$coefficients[,1],
             se= summary(mod)$coefficients[,3])
coef$se <-summary(mod)$coefficients[,4]
coef$UpCI <- coef[,1] + (coef[,2]*2) ### this could be *1.96 but using 2 for simplicity
coef$LowCI <-coef[,1] - (coef[,2]*2) ### this could be *1.96 but using 2 for simplicity

fitted<-data.frame(Open= seq(0,50,2),
               Activity=rep(max(data$Activity),26))

fitted$Marginal <- exp(coef[1,1]*fitted$Open +
                     coef[2,1]*fitted$Open*fitted$Activity)/
              (1+exp(coef[1,1]*fitted$Open +
                     coef[2,1]*fitted$Open*fitted$Activity))

fitted$UpCI <- exp(coef[1,3]*fitted$Open +
                  coef[2,3]*fitted$Open*fitted$Activity)/
             (1+exp(coef[1,3]*fitted$Open +
                  coef[2,3]*fitted$Open*fitted$Activity))

fitted$LowCI <- exp(coef[1,4]*fitted$Open +
                  coef[2,4]*fitted$Open*fitted$Activity)/
            (1+exp(coef[1,4]*fitted$Open +
                  coef[2,4]*fitted$Open*fitted$Activity))