title: What is the difference between sample and outcome? (plus events and observations)

text: I was sure that these are the same things but I do not get the difference reading about mass probability function > Suppose that $X: S → A (A \subseteq R)$ is a discrete random variable > defined on a sample space S. Then the probability mass function $f_X: A → > [0, 1]$ for X is defined as > > $$f_X(x) = \Pr(X = x) = \Pr(\\{s \in S: X(s) = x\\})$$ > > Thinking of probability as mass helps to avoid mistakes since the physical > mass is conserved as is the total probability for all hypothetical outcomes > ''x'': $\sum_{x\in A} f_X(x) = 1$ Well, I just do not get this > The sample space of an experiment is the set of all possible outcomes of > that experiment Sample set, the set of samples is the same as set of outcomes, so we can also call it outcome set. I just do not get why to denote it by two different letters A and S above if these are the same things and why does the function maps the samples to outcomes if it is the same thing? Is outcome a value of the sample? What is the value of the sample then? The functions map values to values. I need a value (of sample) so that function (aka random variable) could produce the value of outcome. * * * **Edit1** I would like to thank who started to elaborate my question. It is not answered, however. I cannot accept it unless the following is resolved. I did not get how events/observations are different from samples/outcomes. Also, I see that you identify samples with _any type_ whereas outcomes are values (you wanted to say _num types_ instead since _any type_ also has values and variables). We, programmers, write functions in terms of arguments/variables. For readability and reliability, we constrain every variable to some type. The values/arguments are placeholders that take concrete values during realization at runtime. The types simply constrain the domain of the variables (types are the domains/set of available values for that variable/argument). I think that it is sorta the same in math, you just skip the name of variable when define function like $f: Int \rightarrow Real$. So, by 'anything', you mean 'value of any, not necessarily numerical, type'. It remains unclear however what is the point of introducing the S domain (aka separating samples and outcomes). You take heads/tails and convert it into a number. Why not to sample real numbers right away (identify S with A, why heads/tails, how introducing S improves our understanding of mass probability function) or why not to make the conversion chain even longer introducing $S_{interm}$ for instance, so that you could sample and use a series of random variables to turn the sample into A, e.g. $S \xrightarrow {randomvar_1} S_{interm} \xrightarrow {randomvar_2} A$? Why not to say that I have 3 pennies and 1 dime and sample A={1,10} right away? Why play around with heads/tails instead? It is also not clear why 'random variables' appear in the S -> A stage rather than in the sampling, to obtain values of S? Does it mean that generation of heads/tails is deterministic whereas mapping them to real domain, {1,10}, is random? Back to the `sample=outcome`. I see that Wikipedia says similarly > A random variable is a real-valued function defined on a set of possible > outcomes, the sample space Ω. So, how can it be that random variable maps samples to outcomes if the outcomes are the domain rather than range of random variables? I think that all this confusion deserves clarification.


title: How to check discriminant analysis assumption in R using lda?

text: I want to use lda in MASS package in R. According to the theory behind that, first need to validate the assumption. Actually I've found some example from the net but they did not bother to validate them. Is lda assumption the same as manova assumption? Am I right? I would appreciate if someone could clarify it.


title: glmnetcr - error in R - regularized continuation ratio logistic regression

text: I'm trying to run regularized ordinal logistic regression with `glmnet.cr()` on 28 predictors, a mix of continuous and categorical. Here are the relevant lines of code:               require(glmnetcr)     attach(data)     newyfac<-as.factor(newyfac)     dim(data)     y<-data$newyfac     x<-data[,2:dim(data)[2]]     mGLMNETCR1<-glmnet.cr(x,y)      `glm()` and `polr()` run fine (with a formula, not "`(x,y)`") but when I try the last line I get this error message:               Error in lognet(x, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  :        NA/NaN/Inf in foreign function call (arg 5)     In addition: Warning message:     In lognet(x, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  :       NAs introduced by coercion      For reference, using `head(data)`, here's what the data looks like:                   newyfac VMONTH  VDAYR   AGE MINOR.  ARRTIME RESIDNCE    SEX ETHUN   RACEUN  ARREMS  PAYTYPER    TEMPF   PULSE   RESPR   BPSYS   BPDIAS  POPCT   ONO2    GCS PAINSCALE   SEEN72  DISCH7DA    PASTVIS RFV1    RFV2    RFV3    EPISODE HOSPCODE     1   4   October Monday  24  FALSE   0.432638889 Private residence   Female  Blank   Blank   No  Medicaid    99.0    64  18  94  61  99  No  Blank   7   No  No  1   Foot and toe pain, ache, soreness, di...    Blank   Blank   Initial visit   1     2   3   October Thursday    15  FALSE   0.854166667 Private residence   Female  Blank   Blank   Blank   All sources of payment are blank    100.0   120 28  138 93  92  No  Blank   0   No  No  1   Labored or difficult breathing (dyspnea)    Blank   Blank   Initial visit   1     3   3   October Thursday    39  FALSE   0.486805556 Private residence   Female  Blank   Blank   Yes Medicare    97.0    98  14  138 88  93  Yes Blank   8   No  No  0   General weakness    Back pain, ache, soreness, discomfort   Disorders of respiratory sound, NEC Initial visit   1      The response, `newyfac`, is ordinal with five levels 1:5.


title: Are randomForest variable importance values comparable across same variables on different dates?

text: Are randomForest variable importance comparable across same variables on different dates? I have a data array X which is of size $T\times N\times K$, where $T=1500$, $N=1500$ and $K=10$. Physically, the 1st index $1,2,\ldots,T$ denotes days, while the 2nd index $1,2,\ldots,N$ represents locations, and the 3rd index $1,2,\ldots,K$ represents the $K$ features/variables measured at each location on each day. The dependent variable is another array $Y$ which is of size $T\times N$. Now I run `randomForest` on each date:               library(randomForest)          importanceValues=matrix(0, T, 10)          for (i in 1:T)          {         y=Y[i, ]              x1=X[i, ,1]         x2=X[i, ,2]         x3=X[i, ,3]         x4=X[i, ,4]         x5=X[i, ,5]         x6=X[i, ,6]         x7=X[i, ,7]         x8=X[i, ,8]         x9=X[i, ,9]         x10=X[i, ,10]              rf=randomForest(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10, importance=T, na.action=na.omit)              importanceValues[i, ]=rf$importance[, 2]          }      As you can see, I can obtain the variable importance values across dates: For example, on the last date above, we have:               > rf$importance                               %IncMSE      IncNodePurity     x1                          311.0803     1113618.9     x2                         4627.7532     3415010.7     x3                         8527.4607     4916842.7     x4                         3507.1872     2919601.3     x5                         2982.0577     2907352.5     x6                         5673.6522     5247811.5     x7                         3893.7793     3618126.4     x8                          135.2311      248212.5     x9                         1759.8080     2334093.9     x10                         852.3294     1562279.1      My questions are:   1. Which one is more useful? (IncMSE or IncNodePurity?)   2. How do I explain to a non-dataminer what "IncNodePurity" is?   3. What's the unit of the "IncNodePurity" column? And can I compare these numbers across dates?    4. One one date, e.g. 9/18/2008, most of the "IncNodePurity" numbers are much larger than those of another date, e.g. 6/1/2012. What can I say about the data sets on these two different dates? (They are different observations for the same variables on different dates) Thank you!


title: Zero-inflated Poisson regression

text: Suppose $ \textbf{Y} = (Y_1, \dots, Y_n)'$ are independent and $$\eqalign{ Y_i = 0 & \text{with probability} \ p_i+(1-p_i)e^{-\lambda_i}\\\ Y_i = k & \text{with probability} \ (1-p_i)e^{-\lambda_i} \lambda_{i}^{k}/k! }$$ Also suppose the parameters $\mathbf{\lambda} = (\lambda_1, \dots, \lambda_n)'$ and $\textbf{p} = (p_1, \dots, p_n)$ satisfy $$\eqalign{ \log(\mathbf{\lambda}) &= \textbf{B} \beta \\\ \text{logit}(\textbf{p}) &= \log(\textbf{p}/(1-\textbf{p})) = \textbf{G} \mathbf{\lambda}. }$$ If the same covariates affect $\mathbf{\lambda}$ and $\textbf{p}$ so that $\textbf{B} = \textbf{G}$, then why does zero inflated Poisson regression require twice as many parameters as Poisson regression?


title: Follow-up updates in case-cohort designs

text: In a case cohort, you randomly sub-sample some feasible number of the cohort (say 20%) and then afterwards add all of the cases that were not part of the sub-sampling. So while a case-cohort is easy to set-up in can be tough to analyse, but there are methods for it. My challenge is how does one deal with follow-up updates in a case-cohort design. For example, if at time A you have 100 cases.. then at time B you have 150 cases. So as more time passes what is a case and what is a control changes. What are the best ways to deal with this? I think one strategy is to redo the sub-sampling, or, just throw in all of the additional cases and sample some new controls with some probability; but this option is not necessarily available. For example, it can be too time consuming to gather all the parameters on the new cases and controls. So the study is no longer a true case-cohort, and more of a nested-case control.. but not really because of the way the controls are sampled. Another obvious answer is to ignore the follow-up, but is there truly no way to harness it? So my challenge is, because of an inability to access the necessary parameters for new cases, if I were to use the follow-up data I would longer have a true case-cohort. But I also am not left with a true nested case-control because my case themselves were not randomly sampled and controls are not matched to the risk level of the cases. Suggestions?


title: When can we add a statistical touch to least square optimization problems?

text: I am trying to connect the dots between statistics and linear algebra/optimization. As you know, Least Square problems are linear algebra and optimization problems. But they also can be connected to statistics. My questions are: Lets say we want to find the $b$'s that minimize the sum of squares of $(y- Xb)$. This is a Least Squares problem. However, we also want to know the "significance" of such $b$'s. It's here that we need to connect the Least Squares problem to statistics. My question is: under what conditions can we do that? We knew that the Gauss-Markov Theorem might play a role here. But the Gauss- Markov Theorem has several conditions. More specifically, if we only care about the "significance" of the $b$'s, what are the "essential" conditions that we have to meet, in order for the $t$-stats to make sense? To be more clear, we don't need the BLUE, we only need to be able to measure the $t$-stats...


title: How does mean and median play an important role in data mining?

text: Can anyone explain me in simple english why does mean and median play a important role on Data mining. Actually finding a mean and median is of what use?   And many people say's median is better at times than mean, why so? And I'm a novice in data mining so it will be really good if the answer is in simple terms. Thanks in advance.


title: What is the distinction between "relationships" and "differences" (particularly in the context of group differences on a numeric attribute)?

text: What is the distinction between a relationship and a difference in the context of statistics - if you have a nominal variable and an interval variable, when should you use a correlation to measure the relationship, rather than a t-test or an ANOVA? I'm not after the usual 'what is your question about - relationships or differences' (although I obviously appreciate the importance), but rather what is it about the data or the constructs that I should know to make this decision. If for example you have males and females and their intelligence scores. Is there a statistical/conceptual reason why a t-test is better here than a point-biserial correlation? This is probably going to be an obvious question. To be honest though, I don't feel like I've ever really understood the distinction, so perhaps some of you can articulate it for me?


title: Using the dependent variable as part of a control variable?

text: I'm regressing on a variable where its ratio with another variable is meaningful in itself. Both are also quite skewed, so I would like to log the ratio so that I have $Y = \beta_0+\beta_1 X_1 ...+\beta\log(\frac{X}{Y})$. It would be nice to have a meaningful coefficient, but unless I'm mistaken it doesn't seem like there is any way to do this. I'm more concerned with whether I can use this as a control variable. Even though it's not strictly linearly dependent with the dependent variable and has another component to be estimated, will this cause problems with the regression?


title: Binary Classification

text: I have conducted 2 experiments on two parameters, $P1$ and $P2$, and I want to classify the two parameters for each experiment. ### Experiment 1               P1:     P2:     0.24    0.14      0.3     0.19      0.37    0.24      0.48    0.29     0.59    0.34      0.67    0.38      0.75    0.43     0.79    0.44      0.8     0.44      0.79    0.45      ### Experiment 2               P1:    P2:       0.13   0.16      0.18   0.2      0.25   0.24      0.31   0.27      0.41   0.32      0.49   0.35      0.56   0.36      0.56   0.35      0.59   0.34      0.54   0.32      What is the easiest way to binary classify these two data sets and what would be the classification accuracy for each experiment?


title: Primitive length measurements

text: A measuring rod has length $u$ (for "unit") and a long object has length $x$. Suppose $u$ is laid end to end $i$ times and $x$ is laid end to end $j$ times. We want to observe whether $iu\ \left\\{\begin{array}{c} < \\\ = \\\ > \end{array}\right\\}\ jx$. But for each iteration of $u$ and of $x$ there is a random error---say we observe whether $$iu + \varepsilon_1+\cdots+\varepsilon_i\ \left\\{\begin{array}{c} < \\\ > \end{array}\right\\}\ jx+\delta_1+\cdots+\delta_j$$ for $i=1,\ldots,I$ and $j=1,\ldots,J$, where the $\varepsilon$s are independent and $\sim N(0,\sigma^2)$ and the $\delta$ are independent of each other and of the $\varepsilon$s and $\sim N(0,\tau^2)$. So we have $IJ$ observations, each binary, equal to either "$<$" or "$>$" (encode them with $0$s and $1$s if you like). **What is known about statistical inference about the ratio $x/u$ in this problem?** Things like the MLE for $x/u$ or the MLEs for $\sigma$ and $\tau$, or confidence intervals for $x/u$, etc. For large values of $I$ and $J$, might one be able to look for things like non-normality of the distributions of $\delta$ and $\varepsilon$?


title: EFA on R: Classifying Observations

text: I am newbie for R and Stats. I am using EFA in R [psych library] to identify the latent factors underlying my data. I got a reasonable amount of factors matching my purpose. My question is once I have the factors, is there a way to categorize observations based on the dominant factors, one or more?


title: Multiple t-tests vs. one-way ANOVA

text: I'm working on a classification problem and I have a very high F1 baseline of 85%. I have trained three classification models and I want to know which one is the best. How can I do so? I tried two ways:   1. To compare each model against the baseline using paired t-test. So I have tests like:      * baseline vs. model 1     * baseline vs. model 2     * baseline vs. model 3 That tells me that only model 1 is significantly higher than the baseline and so I concluded that model 1 is the best. Is this a valid methodology given that usually classification models are compared against baselines?   2. To compare all models in one fell swoop with one-way ANOVA. So I entered the information of models 1-3 and the baseline which gave me a p-value of 0.02 indicating that there is a difference in means. Yet, with pairwise post-hoc tests, there is no significance between any of the pairs. **Which method is the correct one?**


title: Selecting out main effects when locational-scale-invariance isn't an issue(?)

text: I'm starting a new thread to ask a specific question that I'm left with after reading this old, good thread: Do all interactions terms need their individual terms in regression model? The gist of the problem is this: one starts with a complex model motivated by theory. One would like to simplify it -- practically to reduce multicolinearity, and "inferentially" to be able to conduct better inference on the model terms. The standard advice -- well-articulated in the linked thread -- is that it is generally a bad idea to remove main effects of interaction effects, even when AIC improves. This is because the resultant model is no longer invariant to _locational_ shifts in the variables. I emphasize _locational_. But now say you're fitting a model where all the variables on the right hand side have real meanings in the physical world. Lets take, for example, the pH of a substance, or temperature in Kelvin. These things cannot be locationally shifted without changing their meaning. They can be multiplicatively scaled into different systems (i.e.: inches to cm) -- but multiplicative scaling should not affect inference in the ways described in the linked thread. Furthermore, any sort of additive scaling would fundamentally change their meaning -- we'd no longer be talking about the absolute effect of a real thing, but rather the effect relative to some other thing. In my context, I've got a large number of interactions, and I'm agnostic about whether many of the variables should lead to a level effect, or simply moderate the slopes of the responses driven by other variables. For example, I don't know if temperature has its own effect on the phenomenon of interest, or whether it simply affects response to change in pH. So I don't _think_ that I'm making any serious errors in creating an interpretable, plausible model. Would appreciate if anyone could confirm my logic here, or point out any flaws that might be lurking.


title: Is ARIMA better in comparision with Neural Networks?

text: After working on Backpropagation Neural Network and ARIMA Time Series Model, I asked myself which one is better, but can't figure out the answer. They both use different approaches on the same problem (future prediction). Please can someone help me stating the obvious. Backpropagation(C++):                typedef struct {                     /* A LAYER OF A NET:                     */         INT           Units;         /* - number of units in this layer       */         REAL*         Output;        /* - output of ith unit                  */         REAL*         Error;         /* - error term of ith unit              */         REAL**        Weight;        /* - connection weights to ith unit      */         REAL**        WeightSave;    /* - saved weights for stopped training  */         REAL**        dWeight;       /* - last weight deltas for momentum     */     } LAYER;          typedef struct {                     /* A NET:                                */         LAYER**       Layer;         /* - layers of this net                  */         LAYER*        InputLayer;    /* - input layer                         */         LAYER*        OutputLayer;   /* - output layer                        */         REAL          Alpha;         /* - momentum factor                     */         REAL          Eta;           /* - learning rate                       */         REAL          Gain;          /* - gain of sigmoid function            */         REAL          Error;         /* - total net error                     */     } NET;      ARIMA(R):               arima(stockadj,order=c(best.model[1],best.model[2],best.model[3]),xreg=1:n)     stockfor<-predict(stockari,h=100,newxreg=(n+1):(n+100))     ts.plot(stockadj,stockfor$pred,ylab="Original+Predicted Values",main="Forecast")      Thank you.


title: Correct Analysis of Likert score (Product of two Likert scales)?

text: I have two Likert scale questions measuring perception of risk (Severity & Probability). For example: What kind of injury is caused by falling from a height of ten feet?   1. no injuries   2. slight injuries    3. injuries w/medical no permanent damage    4. injury w/medical permanent damage    5. possible fatal injuries  How likely to occur is an accident caused by falling from a height of ten feet?   1. never   2. unlikely   3. somewhat likely   4. likely    5. very likely I am multiplying the responses from these two questions to create a perception of risk score (e.g. Fatal(5) & Never (1) = 5). This computed score lacks a zero (not ratio) and gap between scores are not the same therefore it isn't interval. I consolidate these scores into four groups: acceptable risk or no risk (1 – 4), low risk (5 – 9), moderate risk (10 – 15), high risk (16 – 25) I have successfully conducted Chi-square and Mann-Whitney U tests with the scores converted in the four groupings. But... My question is: Should I be analyzing computed scores themselves? Does the nature of obtaining the score LACK meaning UNTIL I place them in the ordinal 4 categories? Should I ONLY use the consolidated scores (4 groups) in my Chi-square, Mann- Whitney U, and (soon to be done) ordinal regression?


title: Interpretation of probabilities from a mixed-model logistic regression

text: In the following model specification, which is a random intercept 2-level logistic regression: ![enter image description here](http://i.stack.imgur.com/j5NZ0.png) Would two lower level units ($i$) with the same value of $x_{1ij}$ and within the same higher level unit ($j$) but with different values for $y_{ij}$ have the same predicted $\pi_{ij}$?


I have close to a million data sets and whenever I run mean comparison test, either ANOVA or a t-test, I get a significance level of less than .0001 on SPSS. I'm concerned that my sample is so large that of course when I compare the means it will show up as significantly different. Can a sample be too large for ANOVA or a t-test?


title: Chi squared test of homogeneity

text: Suppose you want to test the null hypothesis: H0: p1 = 0.65, p2 = 0.27, p3 = 0.01, p4 = 0.05, p5 = 0.02 to determine whether a population's distribution matches the proposed proportions. You acquire a random sample of 100 individuals. If the null hypothesis were true, what is the expected value of the test statistic? I'm not sure I conceptually/intuitively understand how the chi-squared works. Can someone elaborate on why the answer is 4 here? If the null hypothesis is true, then the population does have proportions equal to 0.65, 0.27, 0.01, etc. So what does this tell us of the test statistic?


title: Confounding variables

text: Suppose we want to study the association between alcohol consumption and cancer. If we include people who also smoke (a confounder), would this lower the relative risk of cancer for alcohol drinkers? If we only include people who do not smoke, this would increase the relative risk of cancer for alcohol drinkers?


title: how to interpret the interaction term in lm formula in R?

text: In R, if I call the `lm()` function in the following way:               lm.1 = lm(response ~ var1 + var2 + var1 * var2)     summary(lm.1)      This gives me a linear model of the response variable with `var1`, `var2` and the interaction between them. However, how exactly do we numerically interpret the interaction term? The documentation says this is the "cross" between `var1` and `var2`, but it didn't give an explanation of what exactly the "cross" is. It would be helpful for me to know what exact numbers R is calculating to incorporate the interaction between the two variables. Thank you very much!


title: Dynamic consistency and multilevel models using lmer

text: I've been using `nlme` and more recently `lmer` to fit multi-level models of time course data using orthogonal polynomials. My colleagues and I originally chose polynomials because we believed that "nonlinear" functions such as the logistic could not be used for multi-level modeling because they are not dynamically consistent. In at least one case this constraint is articulated very explicitly (Willett, 1997, p. 238-239): > In general, the individual growth modeling approach can accommodate any > level-1 model that is _linear in the individual growth parameters_...Many > common growth functions are dynamically consistent, including the quadratic > model cited above and all other polynomial models, regardless of their > order. Other potentially important individual growth models such as the > logistic model (which provides an important theoretical representation of > human development from the perspective of some psychological theories - see > Fischer & Pipp, 1984) is not linear in the individual growth parameters in > its usual formulation. However, I recently discovered that, as I understand it, both `nlme` and `lmer` can use `SSfpl` to fit 4-parameter logistic functions in a multi-level modeling context. Did we misunderstand the dynamic consistency constraint? Perhaps `lmer` and/or `SSfpl` implements the 4-parameter logistic in a dynamically consistent way? If so, does anyone know how it is constrained to be dynamically consistent? Thanks in advance.


title: How to deal with multicollinearity issue when analyzing survey results?

text: I recently conducted my first study looking at parental education as a determinant of risk. I conducted a survey which measured risk in several contexts and also recorded information about the individuals such as gender, height and all the other usual suspects. Little work has been done on the effect of ethnicity on risk and also country of origin/residence on risk (an interesting secondary research question). The issue is that parents nationality are very highly correlated with their children (as one would expect) to the extent that Minitab has discarded various determinants for being highly correlated. Due to the nature of the sample (students in the UK) correlation is something that it unavoidable. I don't wish to discard variables just because they are highly correlated, although I have seen a few suggestions to do so. Is there a way to overcome this issue? Otherwise how else would one look at the marginal effect of country of birth/parents nationality on their risk preferences? Am I overlooking something simple?


title: When is the distribution of product of two normal distributed variables near normal distribution?

text: It is clear the product of normal distributed variables is not normal distributed. For example, if $X \sim N( \mu_1,\sigma_1^2)$, $Y \sim N( \mu_2,\sigma_2^2)$, then $XY$ does not has the distribution of $ N( \mu_1 \mu_2,\mu_1^2 \sigma_1^2+\mu_2^2\sigma_1^2)$. I have been told that even if the distribution of $XY$ is not normal distribution, the distribution of $XY$ is near to normal distribution, when $\mu_1$ and $\mu_2$ are not so small, $\sigma_1$ and $\sigma_2$ are not so big. Is it true? Try following R code:               n1<-rnorm(10000,0,.005)     n2<-rnorm(10000,0,.005)     n<-n1*n2     d<-density(n)     plot(d,lwd=2)     x<-par('usr')     dn<-dnorm(d$x,mean=mean(n),sd=sd(n))     x<-seq(x[1],x[2],length.out=length(dn))     lines(x,dn,col=2,lwd=2)     legend('topright',legend=c('Estimated density','Normal distribution'),lwd=2,lty=c(1,1),col=c(1,2))      ![Density estimation when $\\sigma_1=\\sigma_2=0.005$](http://i.stack.imgur.com/l2R31.png) It seems only when two conditions are both meet, the distribution is near normal. Is there any theoretical analysis?


title: Assigning new observations into existing clusters made from a distance matrix

text: I am building a classification model (binary outcome) and would like to include an external cluster membership code as a predictor. For training the model, this is straight forward. When "scoring" new data with the model, it is less so.   * If the clustering was k-means and used euclidean distance, assigning the new data to its closest cluster (centroid) would be easy.   * **But, how do you do this when you are using a distance matrix as input to the clustering algorithm?** Note: Once the classification model is built, the labeling and meaning of the clusters must remain static. The question is general, but perhaps made more difficult in that I am looking at clustering time series data (using the `dtw` package with `hclust`) - so it seems that one possible solution to the general question ( i.e. to use a predictive model to map new data to the existing cluster assignments) is not possible. Any experience with this?


title: Things to consider about masters programs in statistics

text: It is admission season for graduate schools. I (and many students like me) am now trying to decide which statistics program to pick.   1. What are some things those of you who work with statistics suggest we consider about masters programs in statistics?    2. Are there common pitfalls or mistakes students make (perhaps with regard to school reputation)?    3. For employment, should we look to focus on applied statistics or a mix of applied and theoretical statistics? **Edit:** Here is some additional information about my personal situation: All of the programs I am now considering are in the United States. Some focus on the more applied side and give masters degrees in "applied statistics" while others have more theoretical coursework and grant degrees in "statistics". I'm personally not that intent on working in one industry over another. I have some programming background and know the tech industry a little better than, say, the genomics or bioinformatics industry. However, I'm primarily looking for a career with interesting problems. **Edit** : Tried to make the question more generally applicable.


title: Clustering trajectory patterns with lot of missing values

text: Is there a problem looking for clusters of trajectory patterns in longitudinal data, when much of the longitudinal values were imputed from the baseline data?


title: Delta method and correlated variables

text: I have been reading about the delta method in regards to auto regressive distributed lag models. This is very new to me, so excuse any beginner mistakes. The problem is as follows: We have a model for gasoline consumption. $g$ is the per capita consumption, $y$ disposable income, $p$ is price, $g_{t-1}$ is lagged consumption. All the values are in logs. $$g_t = \alpha_0 + \beta_1 p_t + \beta_2 y_t + \omega g_{t-1} + u_t$$ $\beta_i$ denote the short-run effects and $\beta_i/(1-\omega)$ denote the long-run effect. The problem is that these long-run estimates do not have standard errors calculated in most studies. I found only two papers that do: _Bentzen & Engsted (2001)_ and _Pesaran & Shin (1997)_. They propose to calculate the standard error using the delta method. The problem that I see is that $y_t$ (or $p_t$) and $g_{t-1}$ are highly correlated, thus violating the delta method assumption (as far as I understand it). The correlation is quite clear since both $y_t$ and $p_t$ are significant in the regression above, so taking $$g_{t-1} = \alpha_0 + \beta_1 p_{t-1} + \beta_2 y_{t-1} + \omega g_{t-2} + u_{t-1},$$ we know that there is correlation between $g_{t-1}$ and $p_{t-1}$ (or $y_{t-1})$, given price (or income) persistency, correlation between $g_{t-1}$ and $p_t$ or $y_t$ is surely there. I even dug a whole lot of data from Eurostat to confirm sample correlation and it was there, higher than 0.5 in absolute value. You can also see the estimated standard error using the delta method is much larger than the ones estimated using other methods. That indicates the omitted correlation might cause the overestimation of the standard error. * * * So the question is: Can I use the delta method to estimate the standard error of the non-linear transformation while knowing these variables are correlated? Or does the non-linear nature of the transformation change things? Thanks a bunch.


title: Sample Size for SEM

text: Is there any formula for Structural equation modeling studies?   I know there is a site: http://danielsoper.com/statcalc3/calc.aspx?id=89# Is there any good reference?


title: How to determine if two survey questions are independent?

text: I am trying to determine if there is a dependency between two survey questions about the effectiveness of a subject-specific online communities. The first survey question asks how frequently a user posts questions or comments. The second survey question asks if the participating in the community increases the user's confidence in that subject area. I think I need to to use a Chi-Square Test for Independence. However, I am not sure how to calculate it, especially given the answer choices. The possible answer choices for Question 1 are:   1. Never   2. 1-5 per month   3. 6-10 per month   4. 11+ per month and for Question 2 are:   1. Strongly disagree   2. Disagree   3. Neither agree nor disagree   4. Agree   5. Strongly agree Is the Chi-Square Test for Independence the correct test, and if so, how do I use it against the data? Does the data need to be grouped? For instance, should I group the answers for Question 1 to a binary choice of Never Post or Sometimes Post? Thank you for your assistance and expertise. ==== Further questions.   * I can easily obtain the observed number from the survey data for each cross-section (e.g. Strong Disagree and Never, Agree and 1 to 5 times, etc.). But where do I find my expected number?   * Am I comparing 1 to 5 times, 5 to 10 times, and 11 or more to Never (Never being the baseline or expected amount)?


title: Testing the equality of coeeficients of two different equations in Stata

text: I am working on the joint test using the `Stata`. The test is not of the standard `F-test` form: `y = a*x1 + b*x2` and testing whether `a=b=0`. Instead, the test is of the form: `y1 = a*x` and `y2=b*x` and testing whether `a=b=0.` How do we perform that in `Stata`?


title: How to calculate CDF of g(X)

text: Let $X$ a random variable with distribution $F_X(x)$ $$Y=g(X) = \left\\{ \begin{array}{lr} X-c & : X > c\\\ 0 & : -c < X \le c \\\ X+c & : X \le -c \end{array} \right\\}$$ ## My approach This is a mixed distribution $F_Y(y=0) = P(-c < X < c) = F_X(c)-F_X(-c)$ ## My question: Is this the correct way to write $F_Y(y) $ , $ y \ne 0$. If not, then what is the correct way? $F_Y(y) = P(Y\le y) = P(Y\le y | X>c) P(X>c) + P(Y\le y | X<c) P(X<c)$ This would give. $F_X(y+c) (1-F_X(c)) + F_X(y-c)F_X(c)$


title: Does a time-series have to be stationary before you calculate a z score or t score?

text: It's been a long time since basic statistics. I have a financial time-series that exhibits exponential growth. Before I standardize, do I have to make the time-series stationary? Before I standardize, do I have to normalize?


title: Models for Generalized Estimating Equation?

text: From Wikipedia, Generalized Estimating Equation (GEE) is a method to estimate the parameters of a generalized linear model (with an exponential family distribution for the response). By reading other references online, I am confused whether GEE is an estimation method, or a statistical model like the generalized linear model, but I am inclined to think of GEE as an estimation method, as an alternative to MLE. Am I right? But it only relies on the mean model (as a function of the parameter), and the variance function (as a function of the mean). If I am correct, it is a maximum quasi-likelihood method (although I also want to ask what a maximum quasi-likelihood method is?). So GEE doesn't use the entire likelihood offered by the generalized linear model, but only part of it. Is it correct that GEE should be able to apply to bigger models than the generalized linear model? Thanks and regards!


title: Fitting Weibull regression with ancillary shape parameter in R

text: I am trying to replicate a Weibull model for political leaders' survival in R. The model uses the Weibull distribution with the hazard rate `h(t)=lambda*pt^(p-1)`. The term p is dependent on a variable W. So in `Stata` it can be done simply with adding `anc(W)` to the regression. But the `weibreg` in the `eha` package only allows one to set values to "shape". Is there anyway to do this with `weibreg`? Is there another package/function I should be considering using for this purpose?


title: How to interpret the significant interaction of two non-significant main predictors?

text: I have two non-significant main binary variables in my binary logistic regression model, but their interaction is significant. The variables are centered and no multicollinearity is the case (all VIFs about 1.0). The main variables are nonsignificant, but their interaction is. I want to interpret that significant interaction of two non-significant estimates. It would read something like this: the **effect** of the variable A is less visible in the level 1 of the variable B (B1) and more visible in the second level of B (B2). Or I can say the same thing about the **effect** of B being less visible in A1 and more visible in A2... However, the problem is that neither the effect of A, nor the effect of B are significant! So the above interpretation, although seemingly correct, sounds inconsistent or strange. (How a non-significant effect is supposed to be boosted by the other variable?)... **On a second thought** , it seems that it is actually possible. For example, if I exclude cases with B2 from my sample, now the effect of A in the sample would show up as significant (it is B2 cases in the sample that disallow A to appear as significant)... This is getting more clear in my head now, but do you still have something other than the above interpretation in mind? Thanks a lot.


title: Statistical Significant Difference between 5 groups of data

text: I'm trying to compare 5 different gears used by fishermen and the amount of bycatch found in each. The values have a unit of "% of time encountered" and range between 0% (ie: never have bycatch) and 100% (ie: encounter bycatch every time they go out). I chose this as not all the fishermen fished the same amount of time. I do not have the same number of fishermen using each gear (eg: for gear 1, n = 4 fishermen; gear 2, n=13; gear 3, n=24; gear 4, n= 5; gear 5, n=4). Is there a way of figuring out if there is a statistical significant difference between the bycatch in each gear? I have minitab, but do not know how to use it well yet. Any advice would be very helpful! Thanks!


title: Effect size of the results of an ANCOVA

text: I have compared two regression models using ANCOVA by following this tutorial but using my own data. Here is a snippet of my data:                            X         MC         CC         RC           FM      AT HS     1     0.874375         NA         NA         NA           NA    grab  1     2     0.451250 0.41948802 0.44885230 0.45113473  0.408781437 release  1     3     0.099375         NA         NA         NA           NA    grab  1     4     0.608125 0.50268263 0.52253593 0.51664870           NA release  1     5     0.151875 0.29038872 0.30763473 0.30972255           NA    grab  2     6     0.948750 0.68792116 0.69742615 0.68139072           NA release  2     7     0.452500 0.42250699 0.45583832 0.47252445           NA    grab  2     8     0.894375 0.62101946 0.62549900 0.62338872           NA release  2      MC, CC, RC and FM are four variables that can model X. There are two categorical variables each with two levels: Action Type (AT) and Hand Side (HS). Hand side 1 means right and 2 is left hand. I want to know if hand side is a covariate in any of the models. Consider the case for RC, for this model I perform ANCOVA using the following R command:                aov(X ~ RC * HS, data = X.models)      here are the results:                            Df Sum Sq Mean Sq   F value    Pr(>F)         RC            1 58.427  58.427 10148.368 < 2.2e-16 ***     HS            1 10.732  10.732  1864.025 < 2.2e-16 ***     RC:HS         1  0.069   0.069    11.927 0.0005803 ***     Residuals   859  4.946   0.006         My understanding is that there is significant interaction between HS and RC, and therefore the difference between the slops of the two models is significant. However, in this case the difference between the intercepts is much larger and more interesting. Is there a complementary measurement that would represent the effect size of the interaction? My current solution is to report the difference in slope and intercept alongside the ANCOVA P value, but I am not sure if that is the best way to do it.


title: Determine weights from selection probabilities in a "with replacement" sampling scheme

text: > **Possible Duplicate:**   >  Derive househould weights from a uniformly distributed person sample **EDIT** : Essentially, I have answered the question myself in the linked question: Derive househould weights from a uniformly distributed person sample. I am keeping the text below for reference. * * * For a finite population $U$, I know the probability $\pi_{i,k}$ of sampling case $i$ exactly $k$ times, where $k$ is limited by $m_i$. In my case, $$\pi_{i,k} = \binom{m_i}{k} p^k(1-p)^{m_i - k}$$ for some $p$. I also have a sample of size $n$ (taken _with replacement_ ), but the total size $N$ of the population is unknown. I am trying to derive a reweighting $w_i$ of the sample that will give unbiased and/or minimum variance estimates of totals for, say, an attribute $x_i$. (Estimating $N$ is a special case.) What would be the correct weights for my sample? Is it advisable to remove the duplicates, given that I can reliably identify them? **EDIT** : What if duplicates cannot be reliably identified? My initial guess is $w_i = 1/(1 - \pi_{i, 0})$, but I don't see how to back this with mathematical argumentation. Also, I'd appreciate any hints on literature concerning this topic. (Perhaps this very problem has been treated in a previous paper?) See Derive househould weights from a uniformly distributed person sample for a more practical description of the problem. _The question asked here is different._ Perhaps Calculating % unsampled in sampling with replacement is related.


title: Interaction term using centered variables hierarchical regression analysis? What variables should we center?

text: I'm running a hierarchical regression analysis and I have some little doubts:   1. Do we calculate the interaction term using the centered variables?   2. Do we have to center ALL the continuous variables we have in the dataset, except the dependent variable?   3. When we have to log some variables (because their s.d. is way higher than their mean), do we then center the variable that has just been logged or the initial one? For example: Variable "Turnover" ---> Logged Turnover (because the s.d. is too high compared to the mean) ---> Centered_Turnover? OR would it be directly Turnover --> Centered_Turnover (and we work with this one) THANKS!!


title: Moment Generating Function for Gamma Distribution

text: Let $X$ be a Gamma random variable with shape parameter $\alpha=2$ and scale parameter $\theta=1$. Then the moment generating function of $X$ is $$m_X(t) = \frac{1}{(1-t)^2}, t<1.$$ It is clear that the $t\neq 1$. However, it is also clear that $m_X(t)$ is defined when $t>1$ as shown in the following picture. ![enter image description here](http://i.stack.imgur.com/eUJMA.png) **Then, why do we need to impose the condition $t <1$ here, please?**


title: Confused about GARCH innovation distribution

text: I'm a bit confused about the error distribution of GARCH models. I understand that several standardized distributions (i. e. expected value = 0, variance = standard deviation =1) can be used, for instance standard normal or standardized t. Assuming I use the standardized t distribution and I calculate the standardized residuals, i. e. observations (unless I'm using a nonzero mean term) divided by the estimated conditional standard deviation. Am I right in that I would have to multiply them by $\sqrt{\frac{v}{v-2}}$ where $v$ is the (estimated) df to get something that's (approximately) t-distributed with $v$ degrees of freedom?


title: Testing the influence of monthly fluctuations

text: Let's say we have a dataset that looks like this:               month   year   airquality     jan     2000   36                  ...     ...    ...     dec     2001   33     jan     2002   44     ...     ...    ...     dec     2002   21     jan     2003   44     ...     ...    ...     etc.      The airquality in this case is considered to be the number of parts per million of something. My question is, how would you test if the factor month has an influence on the airquality and how much of an influence? I am having troubles to think of the right method for testing this because besides a monthly difference there may als be an influence of the year factor (ie. airquality may get worse through the years as traffic intensity increases). Hope you can point me in the right direction (the test/model to use, for example)


title: Is there a name for bar charts that replace the bars with color coded objects?

text: I have been trying to find information on this type of chart and just keep getting unintended results. Is there a name or search term to find bar charts that replace the bars with color coded objects like in this example: ![enter image description here](http://i.stack.imgur.com/W1Yer.png) * * * If anyone is interested I'm researching it to offer a solution to this question on the Graphic Design exchange: How do I use designs in a stacked bar chart with precise statistics?


title: statistical analysis for correlation study

text: > **Possible Duplicate:**   >  How to assess association between sleep disorder and a given personality > trait? @Nikita.....Thank You both for responses. I want to do a proposed study on the frequency of a sleep disorder, namely sleep paralysis,and test to see if a trait for mindfulness is correlated with it. A questionnaire will be given to people that may have this disorder to confirm that they do and a scale to measure mindfulness will also be given. So the scale measures for mindfulness are on a 6 point Likert scale so that would be an interval scale I suppose. I guess after I receive the results from the sleep paralysis survey I can categorize the frequencies of the experience accordingly to show the correlation..? So the idea is to simply test the alternate hypothesis that people who experience sleep paralysis more frequently have less trait for mindfulness...So I know this is correlational by design but exactly what statistics to use I am not sure of....I hope this clarifies the intention better


title: Fitting data to the Poisson

text: As a part of my mathematics research, I am trying to model the traffic in my city. I have gathered traffic flow count data from highways and city roads. I am using the Poisson model for distribution. The data I gathered sometimes fails the chi-squared test. In this situation, which other model / addition would you suggest? I want to basically also state the limitations of a particular model in heavy traffic and conclude the best type of model for different traffic flows.   * my data has 13 degrees of freedom ; the chi-square value is coming too large - about 3000 what should i do in this case ?


title: A mixed ANOVA showed no significant interaction, but a dep. t-test showed a significant effect for only one of the groups

text: I have come across something that has got me kind of confused, I hope someone can help me out. I did a mixed ANOVA with one within subjects factor with two levels 'time' and one between subjects factor 'group'. It was a comparison between a positive group and a negative group on scores pretraining and posttraining. The results showed a significant main effect but no significant interaction effect. I should have stopped there maybe, but I did separate dependent t-tests for each group. These indicated that the scores of one group had changed significantly but the scores of the other group had not. So now I don't understand anymore. Wouldn't a significant interaction effect indicate that there was a change in scores and that this change was affected by group (that the decrease in scores was significantly different between the groups)? Since there was no significant interaction shouldn't I find similar changes in both groups (either both significant or both not)? Well, I hope this turns out to be easier than it seems to me right now, merci d'avance. S.


title: Same trend of datasets in R

text: Imagine I have the following data set :               ╔════╦════════════╦═════════════╗     ║    ║ A          ║ B           ║     ╠════╬════════════╬═════════════╣     ║ 1  ║ 0.723185   ║ 0.216422    ║     ║ 2  ║ 0.725254   ║ 0.216104    ║     ║ 3  ║ 0.725166   ║ 0.214883    ║     ║ 4  ║ 0.730094   ║ 0.22244     ║     ║ 5  ║ 0.730336   ║ 0.222342    ║     ║ 6  ║ 0.727938   ║ 0.21453     ║     ║ 7  ║ 0.728334   ║ 0.217435    ║     ║ 8  ║ 0.729766   ║ 0.21452     ║     ║ 9  ║ 0.724968   ║ 0.213434    ║     ║ 10 ║ 0.724352   ║ 0.214565    ║     ║ 11 ║ 0.727014   ║ 0.2212523   ║     ║ 12 ║ 0.731524   ║ 0.22845     ║     ╚════╩════════════╩═════════════╝      How can I test if A and B have the same trend? Is the Jonckheere's trend test really the test I am looking for? I have found this documentation but how would it be in this case?


title: Introductory Text for GAM

text: I'm looking for a text to help someone who uses GLMs in practice become familiar/comfortable with GAMs. Online or physical textbook would be fine. I am approaching this as a practitioner, so I would prefer practical to theoretical, but I'm confident I have the background to digest a more theoretical text if it happens to be particularly worthwhile. I also use R, almost exclusively. If you happen to know of some R oriented texts/tutorials that'd be really appreciated. Thanks, Eric


title: Nonlinear manifold learning and distances between projected points

text: If I use a manifold learning method to project some data points into a low dimensional space, what will be the distances between the projected points? Can I use Euclidean distance? If the distances will not be Euclidean, what kind of distance metric do I have to use and how?


title: Why is svm not so good as decision tree on the same data?

text: I am new to machine learning and try to use scikit-learn(sklearn) to deal with a classification problem. Both **DecisionTree** and **SVM** can train a classifier for this problem. I use `sklearn.ensemble.RandomForestClassifier` and `sklearn.svm.SVC` to fit the same training data(about 500,000 entries with 50 features per entry). The _RandomForestClassifier_ comes out with a classifier in about one minute. The _SVC_ uses more than 24 hours and still keeps running. Why does the **SVC** perform so inefficiently? Is the data set too big for **SVC**? Is **SVC** improper for such problem? Thanks in advanced!


title: Transform hazard rates to other time intervals

text: If a population has a hazard rate of $X$ per year, what would be its hazard rate per month or per day?


title: GLM and GAM package alternative to GRASP in R?

text: I am looking for the best package providing Generalized Linear Models (GLM) and General Additive Models (GAM) for spatial data. The most widely recommended (GRASP) is no longer current in the repository. Is there an alternative that is suggested, or would an archived version of this package still be the best option for these models?


title: Frequency range of wavelet packets decomposition at each level

text: I measured turbulent velocity in 1Hz sampling rate. To remove fluctuations with T>900 sec I decompose the time series by wavelet packets, problem is that I don't know which levels should be kept and which levels should be removed. In other words, if the wavelet packet decomposition is imagined as a tree what would be the frequency range at each node? ![enter image description here](http://i.stack.imgur.com/zp2Ue.png)


title: Strategy for building best fit multiple regression model with time lagged variables

text: I am building a multiple regression model - wrapped in a function - with one dependent variable and a dozen independent variables. The reason why I am building a function is that I need to do this analysis with approximately 75 different datasets. The challenge is that the independent variables correlate better with the dependent variable when they are lagged in time. Unfortunately, not all time lags are the same for each variable and I would like to determine the optimal mix of time lags for each variable while getting the most optimum Adjusted R^2 value for the multiple regression model. Moreover, after building an initial model I will try to reduce the model using the `step(modelbase, direction="both")` function on the model. In the approach I currently have I time lag all the independent variables with the same number of weeks. This results in the best possible model where all independent variables have the same time lag, but I believe (with a valid hypothesis supporting this) that there is a better model out there when we differ the time lag for each independent variable. My question is what is the best strategy to determine the best fit model without making the number of options huge. If I want to determine between 0 and 20 weeks time lag in weekly steps for 12 independent variables I am quickly up to trying to find a match between 4.096e+15 variables (=20^12). I can imagine reducing the problem with the following strategy: Start by finding the best fit model with one independent variable at different time lags. The second step will be to add a second independent variable with its different time lags and find the best model with the two independent variables where the second is tried at different time lags while the first is kept constant. Then add a third variable for which we take a similar approach as the second by keeping the first two variables constant and change try the third with different time lags. Something tells me that this strategy might be decent approach, but something that there also might be a better overall model that contains the not optimal variables for each individual independent variable. Is there anybody who shine some light on how to tackle this challenge?


title: Using first differenced variables in fixed effects model

text: Deal all, It might be a stupid question but since fixed effects estimator either mean centers the data or uses first differences, is it entirely wrong to take first differences of the data and then run fixed effects regression in Stata? Here is the explanation of the problem. So my model is the following $\Delta LP = \alpha+\beta_1 \Delta HC+\beta_2 \Delta FDI+\beta_3 \Delta HC*\Delta FDI$ Where LP is labor productivity, HC is human capital and FDI is foreign direct investment. Since I have panel of 21 countries for 16 years, I need to control for country specific effects and for that Hausman test suggested fixed effects model. However, when we use fixed effects model, it automatically uses first differences of the data. So my question was, is it ok to use first differenced data in fixed effects model? wouldnt it be double first differencing? thank you.


title: Interpretation of Variance and Covariance

text: I am totally new to statistics and I have to create a variance and covariance analysis. I am using SPSS for this. I have created the covariance table: ![enter image description here](http://i.stack.imgur.com/det1L.png) Hopefully it is the right thing. The first 3 variables are the fixed factors, the other 3 are the covariates. Now it is asked from me to determine the "strength of the connection(it probably is a bad translation - I have to determine to what extent there is a connection) of the sum of the square distances of the corrected total variance and of the corrected model". It is a bad translation from German, but I couldn't find an English equivalent and, as I have said before, I do not have any experience with statistics. What is more, how can I determine which variables have a significant influence on the dependent variable?


title: Correlation not significant because there is not enough variance?

text: I have a question about correlations again. I have a dichotomous variable that I want to correlate with the another one (metric) by using the point-biserial correlation coefficient. I get a non- significant result (n=28). But my concern is that one of the variable has almost no variation in the values: 27 observations are $1$, 1 observation is $0$! In other words, no matter what value variable 2 has, variable 1 has in almost every case the value 1. Can this be the reason why a correlation is not significant?


title: Given a kernel, how to find mapping phi?

text: I'm not clear about kernel.   1. How could I construct my own kernel that is valid? Is the only method Merce Theorem (positive semi-definite)? I mean if I know K is a valid kernel, do I know that a*K+b*K, polynomial(K), exp(K), etc., are valid as well?    2. If K is a 2-dimension kernel, how could I expand (create) K to 3-D?    3. Given a kernel, say K(x,y)=(xTy+1)^3, how could I find a phi that make this K = phiT *phi hold?


title: Identification of Integer valued ARIMA process, INARIMA

text: I was wondering if there is any tools like R that can be used to identify the Integer valued ARIMA process? I know that using forecast package in R we can identify the ARIMA process , is there any package for INARIMA? Thank you


title: $\frac{X_T-X_C}{X_T}$ vs $\frac{X_T}{X_C}$ as a metric of treatment effect?

text: I am attempting to understand the use (and misuse) of these two response metrics. When comparing the effect of a treatment to the control:   * one approach is to calculate the relative effect of the treatment:  $$E_1 = \frac{X_T - X_C}{X_C}$$   * another approach is to calculate the ratio of the treatment response to the control response:  $$E_2 = \frac{X_T}{X_C} = E_1 + 1.$$ Do these indices have names? How do their meanings differ? Are there certain situations in which one is preferred to another? * * * _Note: this question has been modified from the original, hopefully providing clarity by removing specific study details._


title: Simulate probit marginal effects with an interaction term

text: I have the following probit model:               library(foreign)     library(MASS)          data <- read.dta("Braith_statecap_JPR.dta")     data <- data[data$year > 1959,]          model1 <- glm(allons3 ~ confbord + rpc + confbord_rpc + neighlgdp + polity2l + polity2sq + lgdp96l + lnpop + postcoldw + peaceall, data = data, family = "binomial"(link = "probit"))      The data along with a stata do file can be downloaded at: Journal of Peace Research: Bratithwaite The paper can be found at: Braithwaite2010 In the paper from which the data is taken, the author simulates the effects of the interaction term in the model (the confbord_rpc term, the multiplication has been done in advance). "confbord" is a dummy variable and "rpc" ranges from 0 to 7. He first calculates the effects with "confbord" set to 0, and "rpc" varying from 0 to 7 (everything else held constant at their means), and then subtracts the effects from the same models with "confbord" held constant at 1. He obtains a very substantial effect from the "cpr" variable through this procedure. I have been replicating this using the following code:                   #simulate coefficients based on the model         beta <- coef(model1)         covvar <- summary(model1)$cov.unscaled         coefs.sim <- mvrnorm(10000, beta, covvar)              #let rpc vary hold confbord constant at 0, and the other variables at their means         constant <- 1         confbord <- 0         rpc <- seq(0,7,1)         ncivwar_rpc <- ncivwar * rpc         neighlgdp <- 8.077         polity2l <- -0.4691         polity2sq <- 56.51         lgdp96l <- 8.123          lnpop <- 9.028         postcoldw  <- 0         peaceall <- 16.14              #bind them together in a matrix            effects0 <- cbind(constant,ncivwar,rpc, ncivwar_rpc, neighlgdp, polity2l, polity2sq, lgdp96l, lnpop, postcoldw, peaceall)              #multiply the effects matrix with the matrix of simulated coefficients to obtain marginal effects         results0 <- coefs.sim %*% t(effects0)              #evaluate each marginal effect on the standardized normal CDF, to obtain predicted probabilities         results0 <- apply(results0,1:2, function(x) pnorm(x))              #do the same procedure for confbord held constant at 1         constant <- 1         confbord <- 1         rcp <- seq(0,7,1)         ncivwar_rpc <- ncivwar * rpc         neighlgdp <- 8.077         polity2l <- -0.4691         polity2sq <- 56.51         lgdp96l <- 8.123          lnpop <- 9.028         postcoldw  <- 0         peaceall <- 16.14              effects1 <- cbind(constant,ncivwar,rpc, ncivwar_rpc, neighlgdp, polity2l, polity2sq, lgdp96l, lnpop, postcoldw, peaceall)         results1 <- coefs.sim %*% t(effects1)         results1 <- apply(results1,1:2,function(x) pnorm(x))              # subtract the estimated probabilities from each other         results <- results1 - results0              #graph the results         x <- seq(0,7,1)         means <- apply(results,2,mean)         sdUpper <- apply(results,2,function(x) quantile(x, .975))         sdLower <- apply(results,2,function(x) quantile(x, .025))         plot(x,means, ylim = c(-1,.2), type = "l")         lines(x,sdUpper, lty = 2)         lines(x,sdLower, lty = 2)         lines(x,rep(0,8), col = "red")      The problem is that with this procedure I get a much lower effect than the author reports in the paper (figure 1), can anyone spot what I am doing wrong?


title: Is it appropriate to compare means using t-test where the raw data are percentages?

text: I want to compare the following percentages:               a<-c(0.59,0.60,0.55,0.55,0.60,0.58,0.67,0.68)     b<-c(0.16,0.21,0.26,0.53,0.84,0.89,0.84,0.94)        * the values for `a` are probabilities calculated with bayes rule   * the values for `b` are number of instances out of the total I did `t.test(a,b)` Is this appropriate? if not how can I compare these two? UPDATE: instead of looking at differences in means is it possible to look at correspondence valu by value? for example compare .59 to .16, then .6 to .21 etc


title: What is the appropriate method to use to calculate customer lifetime value?

text: I'd like to figure out what the potential lifetime value of a customer may be based on their purchasing patterns with our products. I have transactional data that tells me what a customer purchased, when and how much they spent. A customer may have a "lifetime" like so:               Month 1: $300 Purchase License A     Month 3: $400 Purchase License B     Month 12: $150 Renewal of A     Month 15: Do not renew License B     Month 24: Do not renew license A     No licenses left to renew...no longer a customer?          Total Lifetime Value: $850      Customers may purchase as many licenses as they want and the renewal schedule can vary from monthly to yearly. What I'm interested in is finding a proper classification algorithm to differentiate those customers who are not likely to renew versus those who are. For example, if someone renews for two years in a row are they more likely to renew if they've purchased a different product than someone who has only renewed once? I know that for a classification algorithm I need some sort of "output" value to train against. I'm wondering if I can use "total lifetime value" as this output value or some bin of it. Ideally, what I'd like is a method that can take in all these variables and weight them properly. SVGs, Naive Bayes, K-Means all come to mind, but I'm not familiar enough with them to know which one to use or how to properly apply them. If anyone can point me in the right direction I'd really appreciate it.


title: exercise in probability

text: Consider the situation where n items are to be partitioned into $k < n$ distinct subsets. The multinomial coefficients. provide the number of distinct partitions where $n_1$ items are in group $1$, $n_2$ are in group $2$, . . . , $n_k$ are in group $k$. **Prove that the total number of distinct partitions equals $k^n$** I apologise for not using any symbolic notation, don't know how to do it yet. I assume we need to use binomial coefficient to get the left-hand side of the equation i.e $(x+y)^n$. So it basically boils down to proving $x+y = k$.


title: Lower bound for tail of hypergeometric distribution

text: There are several simple and widely used upper bounds on the tail of the hypergeometric distribution, including $P(X > E[X]+tn) <= e^{-2t^{2}n}$, where X is hypergeometric with parameters N, M, and n. (Thinking of the hypergeometric as describing sampling from a population, N is the population size, M is the number of "interesting" items, n is the size of the sample we draw, and X is the number of interesting items in the sample.) This amusing paper is a good summary: > Matthew Skala. Hypergeometric tail inequalities: ending the insanity, 2009. It appears to be an unpublished manuscript, but is available at http://ansuz.sooke.bc.ca/professional/hypergeometric.pdf However, I've been unable to find a simple and reasonably tight lower bound on that same tail probability. Anyone know one?


title: Initialization in pLSA

text: I'm using pLSA (Probabilistic Latent Semantic Analysis). I'm trying to estimate the best value for some parameters with a process of k-fold cross- validation. I have noticed that the model fitting is strongly influenced by the initialization. In order to choose the "best" initialization i'm using the trivial method of multiple random initializations. Does someone know if do it exist a better way?Or maybe some good practices?


title: Latent Dirichlet allocation Implementation

text: I'm looking for some LDA implementation. I know about this one, MALLET but it is coded in Java and I need some more performant. Can someone give me some reference?


title: Prediction intervals predict.Arima r

text: I would like to ask how the long-term (multiple step ahead) prediction intervals are calculated by function `predict.Arima` in R. I am particularly interested in ARIMA models, SARIMA models and in ARIMA models with external regressors (include argument xreg => regression with ARIMA errors)


title: Can you test if one judge gives significantly different scores from other judges?

text: If I had a small number of observations (about 5) in a sample and the population standard deviation was unknown, what test could be done to see if the value of one of the observations was significantly different from the mean? An example of this would be to test if one judge out of 5 judges in a competition was giving a score significantly different from other judges.


title: Controlling for or Excluding in Linear Model

text: I have four variables in my dataset: A volume (numeric) and three factor variables: smoker (yes/no), sports(yes/no) and decease(yes/no). I want to figure out the effect of the factor decease on volume. However, in theory the volumes are supposedly affected by other two factor variables as well. I tried the following model:               lm(volume^(1/3)~decease+smoker+sports)      and the factors smoker and sports were significant while the factor decease was not. However, I have seen some people excluded such observations (smokers and sportsmen), when they are believed to have a predictable effect. I tried that too, having reduced my dataset to about 2/3, and then with the rest of the observations the sole factor decease was significant. Now, which approach should I adopt for the task and why? Thanks.


title: Monotone likelihood ratio property: check my proof; also, who proved it first?

text: I think this is correct: **Lemma**. Suppose $\ell(x,\theta)$ is a (sufficiently regular) likelihood function; then $$\frac{\partial^2\left(\log \ell(x,\theta)\right)}{\partial \theta \partial x}\ge0$$ is equivalent to the monotone likelihood ratio property. **Proof**. (Just symbolic manipulation.) The $x$ subscripts denote the partial derivative w.r.t. $x$. $$\frac{\partial}{\partial \theta}\left(\frac{\ell_x(x,\theta)}{\ell(x,\theta)}\right) \ge0$$ $$\Updownarrow$$ $$\forall x,\forall \theta_1 > \theta_0,\intop_{\theta_0}^{\theta_1}\frac{\partial}{\partial \theta}\left(\frac{\ell_x(x,\theta)}{\ell(x,\theta)}\right) \mathrm{d}\theta \ge0$$ $$\Updownarrow$$ $$\ell_x(x,\theta_1)\ell(x,\theta_0)-\ell_x(x,\theta_0)\ell(x,\theta_1) \ge0$$ $$\Updownarrow$$ $$\frac{\partial}{\partial x}\left(\frac{\ell(x,\theta_1)}{\ell(x,\theta_0)}\right)\ge0 , \mathrm{QED}.$$ Assuming that this is correct, I imagine someone has proved it before. Anyone have a reference?


title: A question on transfer modeling for the intervention analysis of time series data

text: When reading the section of intervention analysis of time series, I have one question regarding the following descriptions. The following graph defines several response patterns for step function $S_t^T$ and pulse function $P_t^T$. ![enter image description here](http://i.stack.imgur.com/0ti2n.png) The step function and pulse function are defined as ![enter image description here](http://i.stack.imgur.com/dEq0B.png) ![enter image description here](http://i.stack.imgur.com/pI0vb.png) With respect to the above patterns, the function of the form (13.1.4) in the following graph is claimed to be equivalent to the one (marked with yellow color in the following graph). Similarly, Equation (13.1.5) is claimed to be equivalent to the other one marked with yellow color. I do not know how to derive this equivalent relationship. ![enter image description here](http://i.stack.imgur.com/Fky7K.png)


title: Differences between logistic regression and perceptrons

text: As I understand, a perceptron/single-layer artificial neural network with a logistic sigmoid activation function is the same model as logistic regression. Both models are given by the equation: $F(x) = \frac{1}{1-e^{-\beta X}}$ The perceptron learning algorithm is online and error-driven, whereas the parameters for logistic regression could be learned using a variety of batch algorithms, including gradient descent and Limited-memory BFGS, or an online algorithm, like stochastic gradient descent. Are there any other differences between logistic regression and a sigmoid perceptron? Should the results of a logistic regressor trained with stochastic gradient descent be expected to be similar to the perceptron?


title: Computing combinations

text: I am having a problem in getting all the combinations from my dataset. I need your suggestions and guidance in doing so. I have a file, whose first column contains the regions and other file five column represents if a binding event takes place or not.               Region  TF1 TF2 TF3 TF4 TF5     Reg1    8   0   4   9   11     Reg2    2   0   0   4   0     Reg3    0   8   9   0   3     Reg4    6   82  8   0   4       Here the most prevalent combination is of TF1,TF3,TF4, TF5 & TF1, TF2, TF3, TF5. . I have a file of 1700 regions (rows) and 25 binding events(columns). How can I get near to my desired result with some associated score for each combination. Thank you.


title: Hypothesis testing terminology surrounding null

text: Which is more correct?   * "Accept the null hypothesis"; or   * "Do not reject the null hypothesis" Is there truly a separate meaning that I'm missing here?


title: General Linear Model (GLM) with Continuous Variable and Categorical Variable (SPSS)

text: I would like to perform a GLM with a continuous variable and categorical variable as fixed factors. For e.g. Weight predicted by Height and Gender. From what it seems, the univariate GLM option in SPSS only allows for categorical variables as fixed factors and not the continuous variable. If I do a linear regression instead, it is the opposite: I can't have the categorical variable as an IV, only continuous. Please help!


title: Count variable as control variable in regression in SPSS

text: I'm doing a research on development of audit fees in 2005-2012. I'd like to see if there's a downward or upward trend in them. I have made a count variable of the years (2005=1 2012=8) and now should use it in the multiple regression to see the trend in the fees. Iknow how to do the multiple regression but what do i so with the count variable in SPSS? Thanks if anyone could help with this!


title: Two-way repeated measures ANOVA with unbalanced data

text: I have a system / operator repeated measures data set for which I would like to answer the questions:   * is system 1 quicker than system 2   * is operator 1 quicker than operator 2   * does the order of using the systems matter So I've hopefully understood the problem correctly to be a "system vs. operator with interactions" question. The measurements in regard to both system and operator are randomised, but unbalanced, i.e., there are different number of measurements for each operator and the order in which they used the systems. (Unfortunate, but the way it goes and somewhat related to the randomisation.) The data is in the format:                   subject time1   time2   order   operator     1   1       1269    422     1       2     2   2       795     327     2       2     3   3       1866    551     1       2     4   4       1263    382     2       2     5   5       1602    438     1       2     6   6       1359    423     2       2     7   7       850     415     2       1     8   8       1080    370     1       2     9   9       568     278     2       1     10  10      582     308     2       1     11  11      1094    654     2       2     ...      I reformatted the data with the make.rm function (http://cran.r-project.org/doc/contrib/Lemon-kickstart/makerm.R) to make it one condition per line. This is so as to allow a "univariate" repeated measures analysis of the data (per the comments of the make.rm function). Therefore the data is now:                   subject order   operator    repdat  contrasts     1   1       1       2           1269    T1     2   2       2       2           795     T1     ...     x   1       1       2           422     T2     y   2       2       2           327     T2     ...      And then I perform the analysis of variance with this data:               analysis <- aov( repdat ~ order * operator + Error( subject / ( order * operator ), data = data.formatted ) )      with a result, but also an error message: "Error() model is singular" So my questions are:   * Have I understood the problem correctly? Am I even doing it right?   * Are the results valid?   * Why do I get an "Error() model is singular"   * What impact does the unbalanced nature of the data play?


title: measurement sampling error and heteroskedasticity

text: I was reading: A potential problem arises from the fact that the measurement (sampling) error associated with the regression beta may vary cross- sectionally. This may introduce heteroskedasticity which would cause inefficient (though unbiased) coefficient estimates. Why would measurement "Sampling" error lead to heteroscedasticity?


title: Why am I getting inconsistent results from Wilcoxon signed rank test and Kruskal-Wallis test?

text: I am doing a study on 3 drugs, comparing response pre-post treatment. My objective is to know if these drugs are effective and which one is better. I used non parametric tests since the results weren't normally distributed and were non transformable. Two drugs were effective with a significant difference on the Wilcoxon signed test and a third showed no significant difference. Yet when comparing the post-results using the Kruskal-Wallis test, no significant difference was observed and the pre-results showed non significant differences between these groups using the Kruskal-Wallis test as well. Why is that? Did I choose the wrong tests?


title: Is NN saturation always bad?

text: I am trying to analyse the effect of hidden unit saturation (outputting mostly 0 and 1 for sigmoid, and not much in-between) on the neural network training performance, and I feel a bit stuck, theory-wise: who says saturation is necessarily a bad thing? The idea is that a saturated unit can not differentiate much between patterns. However, that is only true if the unit outputs the same value for all signals it receives. What if the unit just turns binary? Of course, a collection of binary units has less information capacity than a collection of continuous non-linear functions... But perhaps binary (saturated) units still have a place and should not necessarily be eliminated? What are your thoughts on the dangers of saturation? And... are there any standard cures besides regularisation?


title: Upper/lower standard error makes sense?

text: In the given Image you can see a plot with error bars of different length on each upper and lower side. The bars were generated by calculating the `standard error` once for `all points above` the mean value and once for `those below`. ![http://i.stack.imgur.com/cVlJg.png](http://i.stack.imgur.com/Sq6yW.png) As information on this "method" is very sparse I ask myself wether it is correct to use them. Or is it just the wrong term (lower/upper standard error) I used? And finally what is the meaning of them?


title: Variance of median from the mixture distribution

text: Consider iid samples from a fixed distribution function $F(x)$ and consider its median. Now consider another median from iid samples where one half is drawn from $F_1(x)$ and the other is drawn from $F_2(x)$ with $F_1(x)+F_2(x)=2F(x)$ for all $x$. Which median has a larger variance?


title: How can I clean data of inconsistent format in R?

text: I often deal with messy survey data which requires a lot of cleaning up before any statistics can be done. I used to do this "manually" in Excel, sometimes using Excel formulas, and sometimes checking entries one-by-one. I started doing more and more of these tasks by writing scripts to do them in R, which has been very beneficial (benefits include having a record of what was done, less chance of mistakes, and being able to re-use code if the data set is updated). But there are still some types of data that I have trouble handling efficiently. For example:               > d <- data.frame(subject = c(1,2,3,4,5,6,7,8,9,10,11),     +   hours.per.day = c("1", "2 hours", "2 hr", "2hr", "3 hrs", "1-2", "15 min", "30 mins", "a few hours", "1 hr 30 min", "1 hr/week"))     > d        subject hours.per.day     1        1             1     2        2       2 hours     3        3          2 hr     4        4           2hr     5        5         3 hrs     6        6           1-2     7        7        15 min     8        8       30 mins     9        9   a few hours     10      10   1 hr 30 min     11      11     1 hr/week      `hours.per.day` is meant to be the average number of hours per day spent on a certain activity, but what we have is exactly what the subject wrote. Suppose I make some decisions on what to do with ambiguous responses, and I want the tidied variable `hours.per.day2` as follows.                  subject hours.per.day hours.per.day2     1        1             1      1.0000000     2        2       2 hours      2.0000000     3        3          2 hr      2.0000000     4        4           2hr      2.0000000     5        5         3 hrs      3.0000000     6        6           1-2      1.5000000     7        7        15 min      0.2500000     8        8       30 mins      0.5000000     9        9   a few hours      3.0000000     10      10   1 hr 30 min      1.5000000     11      11     1 hr/week      0.1428571      Assuming that the number of cases is quite large (say 1000) and knowing that the subjects were free to write anything they liked, what is the best way to approach this?


title: Why can't we trust our intuition with probability?

text: If ever there was a case where this become clear is with the Monty Hall problem. Even the great Paul Erdos got fooled by this problem. My question which may be difficult to answer is what is it about probability that we can be so confident of an answer we get uwing an intuitive argument and yet be so wrong. Benford law on first digits and the waiting time paradox are other famous examples like this.


title: Normalization of series

text: I have data compiled by someone else where score averages have been computed over time- averages range from 0-100. The original scores have negative values in many cases and the average would have been negative also, raw average ranges from -30 to 90. How is this 'normalization' accomplished?   Thanks


title: Why (or when) to use the log-mean?

text: I am looking at a scientific paper in which a single measurement is calculated using a logarithmic mean > 'triplicate spots were combined to produce one signal by taking the > logarithmic mean of reliable spots' **Why choose the log-mean?** Are the authors making an assumption about the underlying distribution? That it is what...Log-Normal? Have they just picked something they thought was reasonable... i.e. between a a mean and a geometric mean? Any thoughts?


title: Why is Expectation Maximization algorithm guaranteed to converge to minimum, even local?

text: I have read a couple of explanations of EM algorithm (e.g. from Bishop's Pattern Recognition and Machine Learning and from Roger and Gerolami First Course on Machine Learning). The derivation of EM is ok, I understand it. I also understand why the algorithm coverges to something: at each step we improve the result and the likelihood is bounded by 1.0, so by using a simple fact (if a function increases and is bounded then it converges) we know that the algorithm converges to some solution. However, how do we know it is a local minimum? At each step we are considering only one coordinate (either latent variable or parameters), so we might miss something, like that the local minimum requires moving by both coordinates at once. This I believe is a similar problem to that of general class of hill climbing algorithms, which EM is an instance of. So for a general hill climbing algorithm we have this problem for function f(x, y) = x*y. If we start from (0, 0) point, then only by considering both directions at once we are able to move upwards from 0 value.


title: Dummy / lagged variables and component patterns in regression analysis

text: I've been doing some extra reading for a course I am doing, and I saw this; could someone explain it to me please? _"A time series may exhibit several component patterns. How does regression with dummy variables differ from regression with lagged variables? How do these 2 approaches deal with these component patterns?"_


title: Use deseasonalized sales and combine it with another forecasting method to forecast sales

text: I have sales data that obtain high sales peaks in December, that i deseasonalized such that there is no peaks. I then also applied the Holt- Winter method on the sales. I then combined the deseasonlized sales data with the result from Holt - Winter by taking the average between the two. The peaks as expected give me smaller values then my previous actuals for those periods, but i am not too bothered about it because where i have big dips in the sales due to stock outs, i get good forecasts that fix the gaps in the sales. That is what i am looking for and the peaks will be replaced by the highest value in the history for that point. But i want to know is this method a viable way to do forecasting? Is there something wrong with my method or did someone in the past did something similar. In other words combine decomposed sales data with non-decomposed sales data.


title: Bootstrap of a sample with mixed generative process

text: I want to get importance of a specific predictor of interest in a regression. The data sample is not IID. The sample consists of two sub-samples; twin pairs and singletons. So in order to get a robust estimate of importance I can think of two ways to sample; (1) sample one of each twin-pairs combine it with all the singletons, (2) sample one of each twin-pairs combine it with sample the singletons with replacement. Which way is more appropriate? Or is there any other more appropriate way? I tried both on my data and recorded the F-statistic of anova between the models with and without the predictor of interest (repeated 1000 times). The histogram of (1) looks more normal-like, and that of (2) has a positive skew. A F-test on average F-statistic of course results in (2) giving a higher F-statistic and hence a lower P-value. This might be a naive question but I will appreciate links to statistics literature.


title: Comparing simulated and actual data

text: Say I have a simple linear regression with observed values of y _y = a + bx_ I then have 1000 simulations whereby I generated _y_ randomly. (It's not a matter of reshuffling _y_ ; I generated _y_ through a separate model.) For each simulation I have independent estimates of _b_. My goal is to test whether my actual estimate for _b_ differs from my simulated _b's_. Normally, I would test whether _b_ falls within the 95% confidence interval of the 1000 simulated _b's_. However, each _b_ estimate has it's own standard error. How do I take this into account? I saw one publication that gave the number of simulations that have a _b_ +/i 95% CI that overlaps with the 95% CI of the _b_ from the observed data. Is that the standard way of presenting the results? Also, what if I wanted to show the results (x-y relationship) graphically? What would I use for the 95% confidence interval of the slope of the simulations? Thank you!


title: Calculating optimal number of bins in a histogram for n, where n ranges from 30-350

text: I'm interested in finding as optimal of a method as I can for determining how many bins I should use in a histogram. My data should range from 30 to 350 objects at most, and in particular I'm trying to apply thresholding (like Otsu's method) where "good" objects, which I should have fewer of and should be more spread out, are separated from "bad" objects, which should be more dense in value. A concrete value would have a score of 1-10 for each object. I'd had 5-10 objects with scores 6-10, and 20-25 objects with scores 1-4. I'd like to find a histogram binning pattern that generally allows something like Otsu's method to threshold off the low scoring objects. However, in the implementation of Otsu's I've seen, the bin size was 256, and often I have many fewer data points that 256, which to me suggests that 256 is not a good bin number. With so few data, what approaches should I take to calculating the number of bins to use?


title: Normalizing to zero mean and unit variance before regression

text: I'm new to regression (vector autoregression), and recently encountered the following issue:   1. If I use raw dependent and independent variables to do the regression, the $R^2$, DW-d test and standard error of estimate for each dependent variable are quite good. The $R^2$ is above $95\%$.   2. However, When the dependent and independent variables are normalized to zero mean and unit variables (column normalized), these measures got worse. The $R^2$ ranges from $60\%-90\%$. So what happened after I normalized the data?


title: Fitting a probability distribution to zero inflated data in R

text: I am trying to learn how to fit a probability distribution to a vector of data, using the program R, but there are a lot of potential probability distributions to use! So my question is, how do I find the best distribution for my data, and how do I prove that I have picked the right distribution? Can I acquire AIC values for a whole set of different distributions? The data are observational count data of bees visiting flowers. Each species has a certain number of visits, hence the differing frequencies. The goal is to find the best distribution to describe the bee visitation, show that I have selected the right one, and then use that distribution to sample from randomly for a set of simulations. Here is what the data looks like, it is a vector of count observations. It is zero inflated, with a long tailed distribution (maybe zero-inflated negative binomial?).               i.vec=c(0,63,1,4,1,44,2,2,1,0,1,0,0,0,0,1,0,0,3,0,0,2,0,0,0,0,0,2,0,0,0,0,     0,0,0,0,0,0,0,0,6,1,11,1,1,0,0,0,2)      And here are some basic parameters that I have calculated. I am using standard deviation for sigma, and phi is the proportion of zeroes in the data.               m=mean(i.vec)     #[1] 3.040816     sig=sd(i.vec)     #[1] 10.86078     tab<-table(i.vec)     zero.prop<-as.numeric(tab[1])/sum(as.numeric(tab))     #[1] 0.6122449      As you can see, the standard deviation is much greater than the mean, and I have a very high proportion of zeroes.


title: Calculation of log-logit or log-probit models according to Finney using R

text: I'm trying to implement the **logit/probit model** derivation as introduced by Finney using sample data from http://dge.stanford.edu/SCOPE/SCOPE_12/SCOPE_12.html, chapter 6 (this links to a pdf), page 130, table 6.3. Results are similar but not identical. Finney calculated LD50 = 4.85. Here is the source code and output of my R program:               # Finney, 1952     dosis <- c(2.6,3.8,5.1,7.7,10.2)     nges <- c(50,48,46,49,50)     nok <- c(6,16,24,42,44)          edx <- function(dosis=NA, nges=NA, nok=NA, lk='logit') {       require(MASS)       # weights w according to Finney 1952:       p <- nok/nges       q <- 1-p       w <- nges*p*q            # logit/probit model       edx.data <- data.frame(dosis, nges, nok)       glm.logit <- glm(cbind(nok,nges-nok) ~ dosis, family=binomial(lk), data=edx.data, weights=w)            # Calculation of EDx data       r <- dose.p(glm.logit,p=seq(0.1,0.9,0.2))            # Statistical summary       d <- data.frame(x=c(NA,NA,NA,NA))       rownames(d) <- c('Deg. of freedom','Deviance','1-chi.square','Significant difference between fits and observations')       d$x[1] <- df.residual(glm.logit)           d$x[2] <- deviance(glm.logit)       d$x[3] <- 1-pchisq(d$x[2],d$x[1])           d$x[4] <- ifelse(d$x[3]>0.05,'No','Yes')            # Printing of data tables       print(edx.data)       writeLines('')       print(r)       writeLines('')       print(d)            return(glm.logit)     }          g <- edx(dosis=dosis, nges=nges, nok=nok, lk='logit')      The numerical output is as follows:                 dosis nges nok     1   2.6   50   6     2   3.8   48  16     3   5.1   46  24     4   7.7   49  42     5  10.2   50  44                       Dose         SE     p = 0.1: 1.171955 0.22512512     p = 0.3: 3.615213 0.12064635     p = 0.5: 5.148754 0.09792762     p = 0.7: 6.682294 0.13521132     p = 0.9: 9.125553 0.24565236                                                                                x     Deg. of freedom                                                         3     Deviance                                                 33.3237083095658     1-chi.square                                         2.75205023214653e-07     Significant difference between fits and observations                  Yes      Thus, LD50 would be 5.149 according to this R program (vs. 4.85 as calculated by Finney). What is wrong here? Thank you very much for any help!


title: ML & AI - Algorithm recommendation for "optimization by tuning parameters"?

text: As a newbie in Machine Learning (ML) and Artificial Intelligence (AI), I am here to ask some general algorithm pointers on my following question, which I have simplified into a very generic question. I have a model $M$ defined by a few parameters $p_1, p_2, \dots,p_n$. Now I wish to optimize $M$, i.e., maximize one property value of $M$, by tuning $p_1, p_2, \dots,p_n$. I am looking for algorithms that sort of "learn" what is good and what is bad along the way and find the optimal solution in the end. One similar algorithm in my head is genetic algorithm in AI, but for some reason I cannot use it. So I am here to get some algorithm pointers in ML, hopefully.


title: Forecast model developed now, but only used later, is it still valid?

text: I am not sure what the appropriate title for my question is, so I am open to suggestion/correction. Now, the question:   * I observed a daily series $X_t, t=1,...,T$ also series $Y_t, t=1,...,T$.    * My aim is to use $X$ to predict $Y$ in the future. Specifically, I will observe $X_t, t=T+k,...,T+K$ where $k$ can be large. I will not observe any more $Y_t$. But I would like to estimate it.   * So I fitted a regression with ARIMA error using $(X,Y)_t,t=1,...,T$. (this is done using the `forecast` package.)   * The question is, (a) since new data won't come until some time later ($k$ being large), is any forecast still valid? (b) Even if it is valid, how do I forecast anyway? Do I start from $t=T+1$ and work my away towards $t=T+k$? This seems inefficient because I don't care about $t=T+1,...,T+k$. Furthermore, this isn't really possible because $X$ is not observed during that gap period.(I suppose one can argue you can forecast $X$ using some other model to fill the gap, but this still doesn't solve the problem that it is forecasting things I don't care about, let alone part (a) of my question) So, what should I do? P.S. Extra info. We can assume $T$ is large enough to cover any periods multiple times.


title: Moments of function of Poisson process

text: > Mirror thread on Mathoverflow. (I'm new to Poisson processes, so please edit if my terminology is incorrect.) This is a special case of a problem I'm working on; hoping for intuition that will generalize to a multi-dimensional Poisson process. The context is approximating a distribution defined by the $\arg\min$ of a (non-Poisson) process in Theorem 1 of this paper. **Setup** : consider a Poisson process on $[0,\infty)$, with points $\\{ \Gamma_k\\}_{k=1}^\infty$, $\Gamma_k=E_1+\cdots+E_k$, $E_i\stackrel{iid}{\sim}\rm{Exponential}(1)$. Consider $u=\inf\\{t:1\le\sum_{k=1}^\infty \mathbb 1(\Gamma_k\le t)\\}$, where $1(\cdot)$ is the indicator function (one if true, zero if false). In this special case, $u=\Gamma_1=E_1$, and thus $u\sim\rm{Exp}(1)$. We have a closed form solution for $u$, and the first two (central) moments are $\mathbb E(u)=1$ and $\rm{Var}(u)=1$. **QUESTION** : can we derive $\mathbb E(u)=1$ and $\rm{Var}(u)=1$ _without_ using the closed form of $u$? (In my general problem, there is no closed form.) So we know $u=\inf\\{t:1\le\sum_{k=1}^\infty \mathbb 1(\Gamma_k\le t)\\}$, but we can't just use $u=\Gamma_1$ to calculate moments. **Notes/thoughts** : As a reminder, the mean measure of the considered Poisson process is $m(A)=\lambda(A)$, where $\lambda$ is Lebesgue measure (e.g., $\lambda([a,b])=b-a$; just the total length of set $A\subset [0,\infty)$). This means that the expected number of points in any interval is equal to the length of the interval. Also, the probability of one event occurring in $[t,t+dt]$ is $dt$ since the rate is one here. Also potentially helpful: define the random counting measure $\hat N(t)=\sum_{k=1}^\infty 1(\Gamma_k\le t)$. Then, $\mathbb E(\hat N(t))=\lambda([0,t])=t$, and $u=\inf\\{t:1\le \hat N(t)\\}$. For the first moment, I noticed that in this case, $\inf\\{t:1\le \mathbb E(\hat N(t))\\}=\inf\\{t:1\le t\\}=1=\mathbb E(u)$. That is, solving for $u$ after plugging in the mean measure happens to yield the mean of $u$. This struck me as _not_ true generally: $\mathbb E(u(X))\ne u(\mathbb E(X))$. But maybe some property of the Poisson process and/or $u$'s characterization means that $u$ is a ``linear'' function of the process, so this _is_ generally true here? For the second moment, I haven't made any progress. Tried thinking about the expectation $\mathbb E\left[(u-\mathbb E(u))^2\right]=\mathbb E[(u-1)^2]$, but don't know how I could get that without cheating. Also wondered if there was a ``variance of the process,'' if the plug-in approach for the first moment turns out to be valid. But basically stuck. Last, recall either $\rm{Var}(u)$ or $\mathbb E(u^2)$ is sufficient, since $\rm{Var}(u)=\mathbb E(u^2)-(\mathbb E(u))^2$.


title: Model analysis two-stage least squares

text: I ran a 2SLS model. What does the $\chi^2$ value tell us at the output table? Also, when testing for endogeneity under $H_0$: variables are exogenous   Durbin $p=.25$   Hausman $p=.34$ How can these be interpreted?